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Abstract 

A  new  type  of  hyperspectral  imaging  sensor  is  proposed,  simulated  and  tested, 
which  records  both  spectral  and  2-dimensional  spatial  information.  Dispersive  imag¬ 
ing  spectrometers  typically  measure  multiple  wavelengths  and  a  single  spatial  dimen¬ 
sion.  Unlike  dispersive  imaging  spectrometers,  chromo-tomographic  hyperspectral 
imaging  sensors  (CTHIS)  record  two  spatial  dimensions,  as  well  as  a  spectral  di¬ 
mension,  using  computed  tomography  (CT)  techniques  with  only  a  finite  number  of 
diverse  images.  CTHIS  require  a  reconstruction  algorithm  in  order  to  yield  a  usable 
hyperspectral  data  cube,  and  assume  that  the  point  spread  function  (PSF)  is  known. 
To  date,  the  factors  affecting  resolution  of  these  sensors  have  not  been  examined. 

Lens-based  CTHIS  sensors  use  chromatic  aberration  of  a  lens  and  multiple 
images  in  varying  levels  of  defocus  to  determine  the  chromatic  scene  of  an  object. 
This  type  of  CTHIS  sensor  has  many  practical  advantages  including  simplicity  of  its 
design  and  dual  use  as  a  broad  band  imager  with  no  additional  processing.  The  lens- 
based  CTHIS  concept  has  been  largely  unexplored  up  to  this  time.  The  results  of  this 
research  effort  serve  to  examine  factors  affecting  the  spectral  and  spatial  resolution 
of  a  lens-based  CTHIS  sensor,  specifically  showing  how  many  frames  are  needed  to 
reconstruct  the  spectral  cube  of  a  simple  object  using  a  theoretical  lower  bound.  In 
this  research  a  new  algorithm  is  derived  and  is  used  to  successfully  reconstruct  a 
hyperspectral  object  in  the  presence  of  noise  and  background.  This  new  algorithm  is 
used  to  verify  the  number  of  frames  predicted  from  the  theoretical  bound  calculation 
using  laboratory  data,  thereby  demonstrating  the  validity  of  the  bound  calculation. 
Finally,  a  simple  method  is  proposed  and  tested  to  use  this  sensor  in  the  presence  of 
atmospheric  turbulence  .  This  method  is  shown  in  simulation  to  successfully  remove 
the  effects  of  atmospheric  turbulence  and  estimate  the  atmospheric  seeing  conditions 
blindly  from  raw  lens-based  CTHIS  data. 
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Resolution  study  of  a  Hyperspectral  Sensor 
using  Computed  Tomography 

IN  THE  PRESENCE  OF  NOISE 

I.  Introduction 

1 . 1  Introduction 

Imaging  spectrometers  have  been  used  for  multiple  civilian  and  military  appli¬ 
cations  for  the  last  20  years.  Imagery  obtained  from  these  devices  is  useful  since  it 
contains  both  spectral  and  spatial  information  about  the  scene  under  observation. 
Imaging  spectrometers  measure  the  spectral  content  of  a  light  source,  utilizing  either 
dispersive  optics  (such  as  prism  or  grating  spectrometers)  to  spread  the  spectrum  or 
Fourier  transform  spectra  obtained  by  a  scanning  Michelson  interferometer.  These 
methods  require  scanning  through  either  one  spatial  dimension  (dispersive  spectrom¬ 
eters  -  pushbroom  or  whiskbroom  scanning),  or  the  spectral  dimension  (Fourier  trans¬ 
form  spectrometers)  to  obtain  a  full  spatial-spectral  scene  often  referred  to  as  a  “data 
cube.” 

Recently  a  new  method  for  generating  spectral  imagery  has  been  developed 
which  allows  for  simultaneous  imaging  of  both  spatial  and  spectral  information  using 
computed  tomography  (CT)  algorithms.  Such  imagers  typically  are  called  Computed- 
Tomography  Imaging  Sensors  (CTIS)  or  Chromo-tomographic  Hyperspectral  Imaging 
Sensors  (CTHIS),  although  other  names  have  been  suggested.  CTHIS  use  a  disper¬ 
sive  element  to  project  the  3D  hyperspectral  data  cube  multiple  times  onto  a  single 
2D  image  or  a  few  images  (typically  many  fewer  images  than  required  for  a  pushb¬ 
room  sensor).  Multiple  algorithms  can  be  used  to  take  these  projection  images  and 
reconstruct  the  data  cube  generated  from  other  imaging  hyperspectral  sensors.  Unfor¬ 
tunately,  the  resulting  resolution  of  the  data  cube  can  vary  depending  on  the  system 
setup  and  the  reconstruction  algorithm.  This  research  will  examine  the  effects  of  some 
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different  parameters  on  the  resolution  for  CTHIS,  determine  a  lower  bound  to  predict 
the  number  of  defocus  frames  necessary  to  achieve  a  particular  spectral  resolution  for 
a  lens-based  CTHIS,  propose  a  method  of  reconstructing  CTHIS  in  the  presence  of 
a  large  background,  and  verify  this  reconstruction  method  and  lower  bound  using  a 
laboratory  experiment. 

This  dissertation  is  broken  down  into  six  sections.  The  first  chapter  will  first 
examine  previous  work  in  the  development  of  hyperspectral  sensors,  including  CTHIS 
and  the  basic  overall  system  designs  and  compare  these  sensor  variations.  Chapter 
I  covers  a  description  of  the  background  material  necessary  to  examining  CTHIS 
performance.  Chapter  II  develops  the  theoretical  lower  bound  on  CTHIS  performance 
using  the  Cramer- Rao  inequality  (also  called  the  Cramer- Rao  Lower  Bound  or  CRLB). 
Chapter  III  discusses  the  simulation  parameters  and  setup  of  the  lens-based  CTHIS 
and  develops  a  reconstruction  algorithm.  Chapter  IV  discusses  the  laboratory  setup 
of  the  CRLB  and  how  these  were  matched  to  the  simulation  described  in  chapter 
III.  Chapter  V  discusses  the  results  of  the  CRLB,  the  simulation  and  laboratory 
results  of  the  resolution  measurement  and  how  the  CRLB  can  be  used  as  a  parameter 
predictor  for  resolution  performance  in  the  presence  of  noise.  Finally  Chapter  VI  looks 
at  a  simple  setup  of  a  lens-based  CTHIS  simulation  in  the  presence  of  an  unknown 
atmospheric  Point  Spread  Function  (PSF)  and  a  method  for  blind  estimation  of  the 
atmospheric  point  spread  function. 

1.2  Uses  of  Spectral  Imaging 

Imaging  spectrometers  have  been  used  for  multiple  applications  including  agri¬ 
culture  for  urban  planning,  crop  detection,  mineral  analysis,  chemical  signature  de¬ 
tection,  and  environmental  detection.  The  reflectance  by  various  materials  changes 
with  respect  to  the  wavelength  of  light  incident  on  the  material.  As  light  is  reflected 
from  various  materials,  the  waves  mix  additively  to  form  a  scene  that  may  be  imaged 
by  a  sensor.  Typically  sensors  have  very  broad  bands  for  visual  imagery  grouped  into 
red  (roughly  600  —  700  nm),  green  (500  —  600  nm)  and  blue  (400  —  500  nm).  Hyper- 
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Figure  1.1:  AVIRIS  Hyperspectral  data  cube  of  Moffet  Field 

spectral  imagery  defines  these  bands  more  narrowly  and  in  more  bands  than  typical 
visual  imagery  (as  few  as  5  to  as  many  as  200  bands)  that  are  typically  contiguous. 
Figure  1.1  shows  a  picture  from  the  Airborne  Visual-Infrared  Imaging  Spectrometer 
(AVIRIS)  which  has  approximately  200  bands  from  400nm  to  2.5 fxm  [23].  Figure  1.1 
includes  two  missing  spectral  bands  (the  image  uses  a  pseudo-spectrum  also  known  as 
false-color  imagery).  These  missing  bands  correspond  to  wavelengths  absorbed  by  the 
atmosphere  due  to  water  (top)  and  C02  (bottom).  Also,  note  the  large  red  region  in 
the  upper  right  as  this  corresponds  to  a  high  density  of  shrimp  in  the  pond  near  Moffet 
Field.  Similar  spectral  features  (such  as  the  shrimp)  can  be  used  to  determine  where 
other  materials  of  interest  are.  For  instance,  Figure  1.2  shows  spectral  imagery  col¬ 
lected  from  a  series  of  NASA  satellites  [24]  indicating  the  growth  of  Las  Vegas,  Nevada 
from  a  small  city  in  1975  to  a  much  larger  urban  area  in  2009.  The  light  green  areas 
correspond  to  vegetation  (usually  golf  courses  or  parks),  the  blue  and  dark  gray  areas 
correspond  to  cement  (casinos  and  roads),  the  brown  and  white  areas  correspond  to 
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Figure  1.2:  Spectral  Image  showing  urban  sprawl  in  Las  Vegas,  NV  over  1984-2009 

different  types  of  desert  soil  as  well  as  some  surface  mining  (open-pit)  operations  in 
the  northwest  of  the  city.  This  information  from  these  pictures  can  be  used  by  urban 
planners  to  determine  where  resources  are  or  to  plan  the  next  phases  of  expansion  of 
the  city  as  the  population  grows.  Similar  data  have  been  used  in  rural  areas  to  de¬ 
termine  materials  for  crop  growth  vs  weeds  [29]  and  by  law  enforcement  or  defense  to 
determine  the  uncontrolled  “farming”  of  illicit  drugs  in  open  areas  [15,31].  Figure  1.3 
is  another  example  of  spectral  imagery  (taken  from  [33])  shows  the  east  coast  of  the 
United  States.  Central  Park  can  be  seen  in  the  upper  right  hand  corner  as  a  black  dot 
surrounded  by  green,  which  is  surrounded  further  by  dark  grey  corresponding  to  the 
Manhattan  city  streets  where  the  Hudson  and  East  rivers  come  together  between  Long 
Island  and  Manhattan.  The  black  dot  is  the  Jacqueline  Kennedy  Onassis  Reservoir 
on  the  northern  edge  of  Central  Park.  These  data  also  have  been  used  to  determine 
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Figure  1.3: 


False  color  image  of  Philadelphia,  New  York  and  New  Jersey 
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Figure  1.4:  Combined  Visual/Hyperspectral  image  showing  the  Gulf  of  Mexico  Oil 
spill  2010 

mineral  content  under  the  presence  of  other  materials  (such  as  soil  [4]  or  water  [1]). 
Figure  1.4  shows  an  oil  spill  in  the  Gulf  of  Mexico  in  2010.  The  oil  spill  is  marked 
in  yellow,  while  the  sea  is  shown  with  standard  visual  spectrum  colors  to  show  stark 
contrast  where  oil  is  present.  Figure  1.5  shows  how  the  data  was  collected  to  make  the 
image.  Hyperspectral  data  were  collected  in  20  runs  over  the  Gulf  of  Mexico  using  the 
AVIRIS  sensor.  AVIRIS  can  collect  data  at  wavelengths  from  370 nm  to  2.5 pm  [1]. 
AVIRIS  is  a  type  of  sensor  geometry  called  “pushbroom”  sensing.  The  next  section 
details  how  conventional  hyperspectral  imagers  have  been  designed. 

1.3  Design  of  Hyperspectral  imagers 

Several  designs  for  hyperspectral  imagers  have  been  proposed,  which  fall  mainly 
into  two  categories,  dispersive  imaging  spectrometers  and  Fourier  transform  spectrom¬ 
eters.  Dispersive  imaging  spectrometers  disperse  light  typically  by  using  a  prism  or 
diffraction  grating,  image  a  single  spatial  dimension  and  spread  the  spectrum  across 
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Figure  1.5:  Overplot  of  AVIRIS  flights  with  visual  imagery 
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Figure  1.6:  A  typical  hyperspectral  imaging  sensor 

a  detector  in  the  other  dimension.  Figure  1.6  shows  a  schematic  for  an  imaging  spec¬ 
trometer  using  a  prism.  The  slit  is  imaged  through  the  prism  onto  the  detector.  The 
slit  represents  the  single  spatial  dimension,  and  the  spectrum  is  spread  perpendicular 
to  the  slit  onto  the  detector.  Platform  motion  is  used  to  develop  a  scene  one  line  at 
a  time.  In  a  pushbroom  configuration  (see  figure  1.7  for  a  typical  pushbroom  sensor 
configuration)  as  the  platform  moves  perpendicular  to  the  slit,  the  collection  of  lines 
imaged  through  the  slit  produces  a  hyperspectral  data  cube.  Figure  1.11  gives  a  de¬ 
tailed  sensor  schematic  for  a  pushbroom  sensor.  The  x  dimension  from  Figure  1.7 
is  the  same  as  the  x  dimension  in  Figure  1.11.  A  simulated  hyperspectral  scene  is 
shown  in  Figure  1.9  using  spectra  from  the  Advanced  Spaceborne  Thermal  Emission 
and  Reflection  Radiometer  (ASTER).  ASTER  data  are  available  through  a  NASA 
website  and  can  be  searched  according  to  material  [22],  Figure  1.8  shows  spectra 
for  cement,  glass  and  deciduous  foliage.  An  output  from  a  pushbroom  sensor  for  the 
scene  in  1.9  image  of  line  128  (the  center  of  the  image)  can  be  seen  in  figure  1.10. 
These  spectra  were  chosen  because  they  reflect  large  portions  materials  in  the  image. 
The  glass  in  Figure  1.10  in  the  infrared  range  corresponding  to  8  —  12/im  wavelengths 
although  the  deciduous  foliage  also  mixes  with  the  spectra  in  the  same  area.  De¬ 
ciduous  foliage  has  a  unique  peak  between  3.5  —  5/im  which  can  be  seen  by  yellow 
peaks  in  the  pushbroom  scene  in  Figure  1.10.  And,  hardly  detectable  without  some 
processing  (but  conspicuous  by  its  absence  between  the  windows)  is  the  cement  with 
its  notch  around  lA/im.  The  cement  between  the  windows  can  be  seen  in  the  low 
intensity  areas  around  11/rni  between  the  glass  from  the  windows.  There  are  about 
ten  windows  in  the  scene  across  the  center,  although  the  central  four  are  obscured 
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Figure  1.7:  Pushbroom  Sensing  Geometry 


ASTER  Database  Spectra  of  Basic  Materials 


Figure  1.8:  Simulated  hyperspectral  scene  using  ASTER  spectra 


False  Color  Image  of  Building  with  Foliage 


Figure  1.9:  Simulated  hyperspectral  scene  using  ASTER  spectra 

by  foliage,  and  all  have  metal  window  panes  (seen  by  the  small  lines  in  between  the 
glass  spectra.  Another  method  of  computing  the  spectrum  of  a  scene  involves  using 
interferometers. 

A  Michaelson  interferometer  can  determine  the  spectrum  of  an  object.  Light 
from  a  target  scene  is  transmitted  through  a  beam  splitter  which  sends  the  light 
down  two  equidistant  paths,  with  mirrors  in  each  path  and  then  mixed  resulting  in  an 
interference  pattern  seen  on  the  detector.  A  single  mirror  (called  the  scanning  mirror) 
is  moved  to  produce  constructive  and  destructive  interference.  As  the  interferogram 
changes  along  with  the  motion  of  the  scanning  mirror,  a  profile  of  intensity  is  built  up 
at  each  pixel  on  the  detector.  In  a  Fourier  Transform  Imaging  Spectrometer  (FTIS), 
the  one-dimensional  Fourier  transform  the  intensity  profile  on  each  pixel  is  taken,  and 
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Diffractive  Imaging  Spectrometer  image  of  Single  Slit 


Pixel  Index 

Figure  1.10:  Simulated  hyperspectral  image  taken  by  a  pushbroom  sensor 


the  spectrum  is  derived  by  correlating  this  to  the  position  of  the  scanning  mirror. 
An  example  of  the  intensities  seen  by  a  pixel  is  given  in  Figure  1.13  using  the  same 
line  as  given  in  Figure  1.10.  As  the  mirror  is  moved,  the  intensity  on  each  pixel 
yields  a  spectrum  correlated  to  the  optical  path  length  difference  of  the  two  legs.  The 
center  peak  of  figure  1.13  is  the  maximum  intensity  and  corresponds  to  an  optical 
path  difference  (A opd)  of  0.  The  equation  for  the  intensity  I(Aopd)  seen  by  the  FTIS 
is  given  by: 

poo 

/( Aopd)  =  /  /(A) 

Jo 

The  intensity  per  wavelength  spectra  (A)  seen  by  each  pixel  can  be  reconstructed  by 
using  the  inverse  Fourier  Cosine  Transform  by: 


.  i  2vr 

1  +  cos  [  —A 


opd 


d\ 


(1.1) 


/(A)  =4 


1 \A0pd)  ^I(^opd  0)  COS  ^  ^  A0p^J 


dA 


opd 


(1.2) 


The  dispersive  imaging  spectrometer  has  the  advantage  that  the  spectra  are  imaged 
directly,  and  the  result  can  be  understood  directly  if  there  is  not  a  lot  of  clutter. 
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Figure  1.11: 


Detailed  sensor  schematic  for  pushbroom  hyperspectral  data  cube 
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Also,  the  imager  takes  advantage  of  the  platform  using  the  motion  to  form  an  image, 
which  is  beneficial  for  both  air  platforms  and  satellites.  The  number  of  spectra  is 
fixed,  but  can  be  designed  to  have  relatively  high  spectral  resolution,  however  this 
resolution  is  fixed  with  the  sensor  design  parameters.  The  FTIS  can  be  configured 
to  have  an  arbitrary  spectral  resolution  depending  on  how  finely  the  optical  path 
length  difference  A opci  is  controlled  by  the  sensing  mirror.  This  means  that  the  sensor 
can  be  tuned  for  a  very  high  spectral  resolution,  or  a  coarse  resolution  depending  on 
the  needs  of  the  data  collection.  FTIS  have  a  disadvantage  of  requiring  a  significant 
amount  of  precision  in  the  setup  and  a  stable  platform  to  keep  the  precise  orientation 
of  the  two  paths  of  the  interferometer. 

Another  main  disadvantage  of  these  spectrometers  is  that  they  both  throw  away 
a  significant  amount  of  light.  Dispersive  imaging  spectrometers  disregard  any  light 
outside  of  the  slit  and  therefore  need  a  significant  amount  of  light  in  a  scene.  FTIS 
lose  at  least  50%  of  the  incoming  light,  because  the  light  gets  reflected  out  of  the  front 
of  the  imager.  The  CTH1S  originally  were  designed  to  take  advantage  of  as  much  of 
the  incoming  light  from  a  scene  as  possible. 

1.4  Previous  Work 

1.4-1  Early  Work.  The  use  of  tomographic  imaging  techniques  for  recon¬ 
structing  images  with  two  spatial  and  a  third  spectral  dimension  was  postulated  first 
by  Levin  and  Vishnyakov  [16].  Later,  Okamoto  and  Yamaguchi  [26]  experimentally 
demonstrated  the  first  chromotomographic  sensor  using  a  series  of  amplitude  diffrac¬ 
tion  gratings  yielding  a  diffraction  efficiency  constant  over  all  wavelengths  of  interest. 
Later  that  same  year,  Levin  et  al.  [3]  demonstrated  a  simple  one-dimensional  recon¬ 
struction  using  prisms  with  variable  dispersion  believing  they  were  the  first  to  demon¬ 
strate  the  concept.  Okamoto,  et  al.  [25]  used  phase  diffraction  grating,  yielding  even 
more  light,  at  the  expense  of  a  more  complex  reconstruction  technique  because,  in 
this  case,  diffraction  efficiency  is  dependent  on  wavelength.  Compared  to  techniques 
developed  by  later  work,  simple  algebraic  reconstruction  techniques  were  used  similar 
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FTIS  image  seen  by  detector 


Figure  1.13:  Detailed  sensor  schematic  for  pushbroom  hyperspectral  data  cube 

to  other  computed  tomography  problems.  However,  both  groups  demonstrated  the 
advantages  of  a  chromotomographic  sensor’s  potential  for  high  optical  throughput 
and  the  ability  to  image  flash  events  over  other  spectral  sensors.  The  disadvantage 
of  these  designs  is  that,  while  they  are  simple,  they  do  not  offer  significant  spectral 
resolution  when  compared  with  dispersive  imaging  spectrometers  or  FTIS. 

1-4-2  Crossed  Phase  Gratings.  Descour  and  Dereniak  [7]  detail  experimen¬ 
tal  results  of  a  CTHIS  using  a  series  of  crossed  phase  gratings.  Their  work  can  be  seen 
as  extended  the  work  of  Okomoto  and  Yamaguchi  to  an  even  larger  number  of  wave¬ 
lengths  by  using  a  statistical  reconstruction  technique.  The  authors  further  developed 
the  theory  of  CT  sensing  using  a  discrete-to-discrete  method  for  reconstructing  the 
object  cube  that  takes  into  account  noise  sources  and  diffraction  efficiency.  In  [7],  a 
filter  was  used  to  reduce  the  spectral  range  of  the  orders  allowing  multiple  diffraction 
orders  to  be  detected,  because  these  multiple  orders  correspond  to  differing  angle  pro¬ 
jections  of  the  hyperspectral  data  cube  and  therefore  contain  different  information. 
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In  [8],  the  authors  use  a  computer  generated  hologram  as  a  dispersion  element  which 
is  specifically  designed  for  diffraction  efficiency  and  image  location.  The  advantage  of 
using  diffraction  gratings  and  holograms  specifically  designed  for  CTHIS  extended  the 
possibilities  of  this  sensor  to  real-world  applications  rather  than  a  simple  laboratory 
setup  [14],  and  also  that  these  sensors  can  image  flash  (one-time)  events.  One  of  the 
disadvantages  of  using  diffraction  gratings  is  the  inconsistent  spectral  efficiency  at  all 
wavelengths  (the  sensor  will  be  more  sensitive  at  some  wavelengths  than  others).  This 
was  the  first  demonstration  in  the  literature  of  a  real-world  use  of  a  CTHIS,  however 
as  referenced  in  [14]  there  were  significant  issues  that  needed  computer  processing  to 
recover  the  actual  imagery,  and  some  of  this  was  directly  due  to  the  design  of  the 
sensor. 

1-4-3  Direct  Vision  Prism.  In  [21],  Mooney  details  a  CTHIS  design  differing 
from  earlier  work  using  a  rotated  direct  vision  prism  in  order  to  sample  the  spatial- 
spectral  object  cube.  All  earlier  work  used  a  diffraction  grating,  or  a  series  of  prisms 
with  different  dispersions  in  order  to  change  the  angle  sampling  the  object  cube.  The 
DV  prism  differs  from  a  standard  prism  in  that  the  central  wavelength  passes  through 
without  any  angular  deviation  [20,21],  while  a  known  angular  dispersion  is  applied 
to  other  wavelengths.  Using  a  DV  prism  also  opens  up  an  interesting  opportunity  to 
calculate  vector-based  images,  similar  to  a  linear  dispersive  imaging  spectrometer  [12], 
This  method  differs  from  a  standard  CTHIS  (using  a  DV  prism)  only  by  software  and 
can  be  useful  in  quickly  determining  scenes  of  interest  which  can  then  be  imaged  using 
the  standard  DV  prism  CTHIS  reconstruction.  The  disadvantage  of  these  sensors  is 
that  they  require  multiple  snapshots  of  a  scene  to  be  effective  whereas  the  diffraction- 
based  CTHIS  sensors  can  image  flash  events. 

1-4-4  Chromatic  Lens  Aberration.  In  [17]  Lyons  proposed  the  use  of  a 
diffractive  optic  element  specifically  designed  to  focus  wavelengths  to  varying  distances 
similar  to  a  Fresnel  lens.  This  research  used  simple  images  directly  to  determine  the 
spectral  content  correlating  that  with  the  position  of  the  lens.  This  element  yields 
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a  small  change  in  the  focus  distance  yielding  large  change  in  the  focal  length,  but 
has  the  disadvantage  of  having  a  large  changes  in  diffractive  efficiency  at  different 
wavelengths.  In  [5]  Cain  also  proposed  using  the  chromatic  aberration  of  a  lens 
as  a  dispersion  element  and  moving  the  images  in  and  out  of  focus  for  successive 
wavelengths  to  capture  multiple  images.  This  has  the  advantage  of  a  simple  optical 
setup  rather  than  using  expensive  optics  for  the  dispersion  element.  Also,  using  this 
method,  an  existing  optical  system  can  be  turned  into  a  CTHIS  using  only  slight 
modifications,  such  as  the  addition  of  a  telescoping  lens  and  an  aperture  stop.  The 
disadvantage  to  this  method  is  that  the  magnification  changes  with  respect  to  each 
wavelength  may  need  to  be  accounted  for  in  the  reconstruction  depending  on  the  lens. 
Also,  this  design  suffers  from  the  same  disadvantage  of  direct  vision  prisms  requiring 
multiple  frames  in  order  to  compute  a  spectral  scene. 

1.5  CTHIS  sensor  design 

The  basic  motivation  for  CTHIS  is  twofold:  first,  to  make  use  of  all  available 
light  and  second,  to  reduce  or  eliminate  the  amount  of  scanning  necessary.  Typical 
dispersive  imaging  spectrometers  make  use  of  a  slit  which  is  re-imaged  onto  the  detec¬ 
tor.  Figure  1.6  shows  the  design  of  a  standard  dispersive  imaging  spectrometer.  This 
slit  is  typically  the  size  of  a  single  line  of  pixels.  A  typical  CTHIS  detector  replaces 
the  slit  with  a  wider  held  of  view  held  stop  (figure  1.14),  and  the  dispersion  element 
with  either  a  specially  made  diffraction  grating  or  a  direct  vision  (DV)  prism. 

Two  methods  can  be  used  for  capturing  spatial-spectral  frames  for  computed  to¬ 
mography,  single-frame  and  multi-frame  detection.  Single-frame  designs  for  a  CTHIS 
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Figure  1.15:  (a)  a  prism  projects  an  object  along  a  single  prism  axis,  (b)  a  grating 

projects  an  object  along  multiple  axis  simultaneously 

differ  significantly  from  dispersive  imaging  spectrometers  in  that  the  gratings  are  de¬ 
signed  to  diffract  light  in  two  dimensions  and  over  multiple  diffraction  orders.  A 
single-frame  design  also  has  the  added  advantage  of  being  able  to  capture  all  wave¬ 
lengths  at  a  single  time.  In  [7],  a  filter  was  used  to  reduce  the  spectral  range  of 
the  orders,  allowing  multiple  diffraction  orders  to  be  detected,  because  these  multiple 
orders  correspond  to  differing  angle  projections  of  the  hyperspectral  data  cube  and 
therefore  contain  different  information.  The  gratings  (or  other  suitable  diffraction  el¬ 
ement)  designed  for  each  of  these  configurations  [7,8,34]  were  chosen  to  maximize  the 
amount  of  light  in  each  diffracted  order  used  and  to  maximize  the  detector  area  used 
to  capture  CTHIS  frames.  However,  since  the  goal  is  to  spread  multiple  orders  over 
the  entire  detector,  diffraction  efficiency  is  not  constant  with  respect  to  wavelength 
or  order. 

In  multi-frame  detection,  a  DV  prism  is  designed  such  that  some  center  wave¬ 
length  Ao  (or  the  whole  multi-color  scene)  passes  through  without  deviation.  Multiple 
frames  are  collected  by  rotating  the  prism  along  the  axis  of  the  undeviated  wavelength 
(figure  1.15(a)  [21]).  DV  prisms  are  generally  manufactured  from  2  prisms  with  dif- 
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ferent  indices  of  refraction  which  cancel  out  the  angular  dispersion  at  the  center 
wavelength.  The  object  and  viewing  conditions  are  assumed  to  be  constant  during 
the  rotation  of  the  prism.  The  number  of  frames  required  can  be  much  fewer  than 
those  required  for  scanning  of  other  spectrometers.  Pushbroom  dispersive  imaging 
spectrometers  require  the  number  of  frames  to  be  the  spatial  pixels  in  the  scanning 
direction  and  Fourier  transform  spectrometers  require  the  number  of  frames  to  vary 
as  the  total  bandwidth  of  interest.  Also  of  note  is  that  keeping  a  single  prism  aligned 
should  be  much  simpler  than  keeping  two  arms  of  an  interferometer  aligned,  so  the 
scanning  of  a  DV  prism  CTHIS  is  a  reduction  in  complexity  over  a  Michelson-based 
Fourier  transform  spectrometer.  CTHIS  scanning  is  only  required  to  adequately  sam¬ 
ple  the  projections  of  the  data  cube.  Figure  1.16  shows  an  example  of  a  DV  prism 
as  the  rotation  varies.  Note  that  in  this  picture  A2  represents  the  undeviated  wave¬ 
length.  As  will  be  seen  in  the  algorithms  discussion,  there  is  inherently  some  missing 
information  in  CTHIS  data  which  must  be  reconstructed.  This  missing  information 
can  be  traded  off  with  processing  time  and  estimated  vs.  measured  resolution.  The 
main  disadvantage  of  using  a  DV  prism  is  that  they  are  difficult  to  manufacture,  and 
may  be  expensive  depending  on  the  wavelengths  of  interest. 

Another  method  using  a  multi-frame  design  is  to  use  the  chromatic  aberration 
of  a  lens  to  achieve  dispersion.  Due  to  chromatic  aberration,  the  focal  length  of  the 
lens  will  be  dependent  on  wavelength.  A  central  wavelength  for  the  lens  is  chosen, 
and  the  lens  position  is  varied  around  the  corresponding  focal  length.  Note  that  in 
figure  1.15  that  the  dark  square  represents  a  single  undeviated  wavelength  for  the 
prism  in  (a),  and  corresponds  to  a  sum  over  all  wavelengths  using  the  grating  in  (b). 

1.6  Reconstruction  Algorithms 

Having  discussed  the  basic  design  of  a  CTHIS,  we  now  turn  our  examination 
to  the  reconstruction  of  the  data  cube.  The  processing  is  an  essential  step  to  using 
CTHIS  imagery  because  without  it,  the  spatial  and  spectral  data  are  multiplexed  in 
the  image.  The  reconstruction  mathematics  borrow  heavily  from  the  already  estab- 
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Figure  1.16:  Example  of  a  DV  prism  in  a  CTHIS  sensor 

lished  world  of  computed  tomography,  using  many  of  the  mathematical  formulations 
applied  to  medical  imaging,  Synthetic  Aperture  Radar  and  other  applications.  Al¬ 
though  some  of  the  mathematical  underpinnings  from  may  be  borrowed  from  CT, 
some  parts  of  the  reconstruction  problem  still  require  special  considerations  in  order 
to  accurately  reconstruct  the  data  cube.  CTHIS  suffer  from  a  fundamental  limitation 
in  that  the  projection  into  the  image  plane  cannot  fully  cover  the  information  neces¬ 
sary  to  directly  reconstruct  the  data  cube.  Methods  have  been  proposed  to  overcome 
this  limitation  based  on  principle  component  analysis  (PCA)  and  projections  onto 
convex  sets  [20,21],  We  will  first  look  at  the  reconstruction  problem,  followed  by  the 
problem  areas  associated  with  reconstruction.  Then  we  will  discuss  reconstruction  al¬ 
gorithm  types,  and  finish  this  section  with  a  discussion  of  comparisons  and  contrasts 
between  the  different  algorithms. 

1.6.1  Algebraic  Reconstruction.  In  [25],  Okamoto  and  Yamaguchi  develop 
the  equations  for  the  images  seen  by  a  single  frame  grating-based  CTHIS.  Using  x,y 
as  index  variables  in  the  imaging  plane,  A  as  the  wavelength,  o(u,v,  A)  as  the  object, 
and  u,v  as  spatial  coordinates  in  the  object  plane,  an  image  seen  by  the  detector  is 


19 


given  by: 


i(x,y,  A) 


o(u,  v,  A )h(x  —  u,y  —  v,  \)dudv 


(13) 


J  u  J  V 

where  h(x,  y ,  A)  is  the  point  spread  function  of  the  imaging  system.  In  this  case,  the 
standard  grating  equation  can  be  used,  and  yields: 


h(x,  y,  A) 


Y  VjkWS  (x-  —  k,y- 

j,k=- 1,0,1  '  1 


(1.4) 


where  rj3  j.  are  the  diffraction  efficiencies  of  the  jth  and  kth  diffraction  orders,  A  is  the 
grating  constant,  and  5  is  a  Dirac  delta  function.  In  [26],  only  the  —1,  0, 1  diffraction 
orders  were  considered.  Finally  the  noise- free  signal  detected  by  the  imager  is  given 
as  a  simple  example  by: 


i(x,y)  —  f  a(\)i(x,y,  X)d\  —  Y  %fc(A)°  (x  ~  V  ~  (L5) 

J  j,k=- 1,0,1  A  ' 

where  rj'jk  is  the  product  of  diffraction  efficiency  rjjk  and  corresponding  detector  spec¬ 
tral  response  a(A).  The  reconstruction  given  in  [25]  is  a  modified  algebraic  recon¬ 
struction  technique  (MART)  which  is  an  iterative  technique  useful  for  small  number 
of  projections.  Setting  an  initial  object  estimate  o°(x,y,  A)  to  a  positive  constant, 
successive  object  estimates  are  then  calculated  by  iterating  over: 

og+\x,y,  A)  =  °q(x,y,  A)  (1.6) 

where  iq(x,y )  is  the  projection  of  the  <yth  estimate  of  the  signal  and  is  recalculated 
using  equation  (1.5)  with  the  gth  estimate  of  the  object.  This  algorithm  is  similar  to 
the  maximum  likelihood  estimator  given  in  section  1.6.3. 


1.6.2  Projections  onto  Constraint  Sets.  If  the  optical  setup  is  based  on  the 
DV  prism  as  in  [20,21],  the  formulation  of  the  reconstruction  becomes  a  convolution 
of  multiple  wavelengths  with  respect  to  each  rotated  frame.  In  [21],  a  DV  prism  setup 
with  three  separate  wavelengths  and  four  separate  frames  the  forward  model  is  given 
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by: 


h(x,y)  =  hltl{x,y)  *  ox{x,y)  +  hh2(x,y)  *  o2(x,y)  +  h1>3{x,y)  *  o3{x,y) 

12 (x,y)  =  h2,i(x,y)  *  oi(x,y)  +  h2,2(x,y)  *  o2(x,y)  +  h2,3(x,y)  *  o3(x,y) 

(1.7) 

k(x,y)  =  h3>1(x,y)  *  oi(x,y)  +  h3,2(x,y)  *  o2(x,y)  +  h3}3(x,y)  *  o3(x,y) 

U (x,y)  =  hAil(x,y)  *  Oi(x,y)  +  h±2(x,y)  *  o2{x,y)  +  h^3(x,y)  *  o3(x,y) 

Where  om(x,y )  is  the  spectral  distribution  corresponding  to  band  m,  ik(x,y )  is  the 
data  recorded  for  the  prism  orientation  k,  hk,m(x,y )  is  the  point  spread  function  for 
spectral  band  m,  with  respect  to  prism  orientation  k  and  the  convolution  operator  is 
represented  by  the  *.  The  summation  form  of  (1.8)  is  given  by: 

N 

ik{x,  y)  =  hk,m{x,  y)  *  0m(x,  y)  (1.8) 

m=  1 

A  geometrical  optics  model  for  hktm(x,y )  [21]  is  given  by: 

hk,m{x,  y)  =  5[x  -  (k  -  k0 A  cos {<j>m)),y  -  (k  -  k0A  cos (</>m))]  (1.9) 

where  is  the  index  of  the  initial  prism  orientation,  A  is  an  offset  defined  by  the 
prism  geometry,  and  is  the  angle  for  spectral  band  m  defined  by  the  prisms  relative 
index  of  refraction.  In  the  Fourier  domain,  the  rotation’s  spatial  spectra  are  a  product 
of  the  original  spatial  spectra  multiplied  by  the  optical  function.  Taking  the  Fourier 
transform  representation  of  (1.8)  using  matrix  notation  becomes: 

I(Z,0  =  H(£,QO(£,0  (1-10) 

where  the  O,  /,  and  H  are  the  Fourier  transforms  of  o,  i,  and  h  respectively.  The 
problem  then  becomes  a  larger  number  of  inversions  of  smaller  matrices  which  is  the 
motivation  for  separating  the  projections  over  multiple  prism  rotations  rather  than  a 
larger  matrix  inversion  necessary  for  a  single-frame  approach.  For  instance,  in  [21], 
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o(x,y,A) 

y 


i 


Figure  1.17:  Example  of  3-D  object  reconstruction  using  Fourier  transforms 


the  authors  solve  “50  240x240  spectral  images  from  80  240x240  images  using  57600 
80x50  matrix  inversions  rather  than  a  single  4,608,000x2,880,000  matrix  inversion” 
necessary  for  the  same  problem  using  a  grating  rather  than  a  DV  prism.  This  formu¬ 
lation  directly  correlates  the  spectral  Fourier  transform  with  the  measured  data.  The 
Fourier  transform  exhibits  a  special  property  known  as  the  projection-slice  theorem, 
the  Fourier  slice  theorem  or  the  central  slice  theorem,  all  of  which  state  that  the 
Fourier  transform  of  a  projection  can  be  used  as  slices  through  a  higher- dimensional 
Fourier  space  (see  figure  1.17). 

The  H  matrix  and  its  pseudoinverse  can  be  precomputed  and  stored.  The 
pseudoinverse  of  H  can  also  be  diagonalized  using  the  singular  value  decomposition 
and  the  forward  and  inverse  Fourier  transform  can  be  directly  applied.  Inverting  the 
H  matrix  is  less  complex  computationally  than  inverting  the  similar  matrix  mentioned 
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in  [7]  because  the  matrices  can  be  diagonalized  due  to  the  layout  of  the  DV  prism 
CTHIS.  Diagonalizing  H  yields: 


H  =  UWV] 


(1.11) 


where  W  is  a  weighted  diagonal  matrix  according  to  equation  1.12,  and  U  and  V 
are  the  weighting  matrices  needed  for  the  product  to  equal  H.  This  is  also  the 
singular  value  decomposition  (SVD)  of  H ,  where  U  and  V  are  unity  when  multiplied 
by  their  transpose  conjugate  (also  called  the  Hermetian  adjoint).  Inverting  (1.11)  and 
substituting  back  in  to  equation  (1.10)  we  find  equation  (1.13)  in  which  the  inverse 
of  W  (elements  referred  to  as  w^j))  is  replaced  by  W~l  (elements  referred  to  as  whj) 
where  the  small  singular  values  elements  of  W~x  are  replaced  with: 


WiAt’Q 


«&(&  0  +  v2 


(l.!2) 


so  as  not  to  amplify  noise.  In  this  case,  rj  is  a  tuning  parameter  (dependent  on  noise, 
data  and  sampling)  generally  close  to  unity. 


O  =  VW~lUU  (1.13) 

In  practice  however,  the  direct  matrix  inversion  (or  pseudo-inverse)  is  not  solely 
the  solution  to  the  reconstruction  problem,  but  it  is  used  as  an  estimate  to  determine 
the  input  to  an  iterative  algorithm  which  then  calculates  the  reconstruction.  We  can 
make  use  of  the  relationship  between  the  Fourier  domain  and  the  captured  images  [2] 
to  construct  a  simple  iterative  algorithm.  For  the  rest  of  this  development,  we  assume 
that  we  can  reorganize  the  two-dimensional  matrices  as  single-dimensional  vectors  to 
simplify  notation.  By  adding  the  forward  and  inverse  Fourier  transforms  to  H,  adding 
noise,  and  expanding  the  psuedoinverse,  equation  (1.8)  becomes: 

i  =  F~1UWV*Fo  +  n  (1.14) 
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where  n  is  a  noise  vector,  and  F,  F  1  are  the  forward  and  inverse  Fourier  transforms 
respectively.  Our  iterative  estimate  of  o  becomes: 


Oj  =  F-'FlLV-'rfFi  +  Rv*FFoj-i 


where  L  and  R  are  diagonal  with  elements: 


Uj  «,o 


WjK.O 
Wy  (f,  0  +  P2 


KyK.C)  =  !-£««,<) 


(1.15) 


(1.16) 

(1.17) 


This  algorithm  doesn’t  directly  take  into  account  fixed  pattern  noise  of  the  focal 
plane  array  (FPA),  however  the  authors  note  that  typically  these  errors  get  mapped 
to  the  undeviated  (central)  spectral  band.  In  [20]  Brodzick  and  Mooney  further 
generalize  the  model  developed  in  [21]  and  note  that  it  is  of  a  class  of  algorithms 
called  projections  onto  convex  sets  (POCS).  One  POCS  algorithm  that  is  also  well- 
known  is  the  Gercheberg-Papoulis  algorithm,  and  is  similar  to  the  generalization  of 
the  algorithm  mentioned  in  [21].  Their  derivation  amounts  to  projecting  the  Fourier 
domain  along  the  principle  component  axes  of  the  range  space  (the  known  portion) 
of  the  system  transfer  function  matrix  H  (called  A  in  [20]  and  P  in  [21])  and  the 
projection  of  the  unknown  portion  into  the  null  space  of  H .  This  null  space  projection 
is  called  the  transform  domain  constraint.  Then  the  algorithm  uses  the  fact  that  the 
object  domain  (the  hyperspectral  data  cube)  has  a  high  degree  of  correlation  between 
bands.  This  redundancy  can  give  an  advantage  by  using  the  PCA  projection  of  the 
object  domain  along  the  principle  components  to  fill  in  the  unknown  portion  of  the 
frequency  domain.  By  using  PCA,  the  algorithm  is  forcing  the  Fourier  domain  data 
to  be  consistent  with  the  object  domain  principle  component  projection.  This  second 
constraint  is  called  the  object  domain  constraint  and  is  a  generalization  of  [21].  Small 
singular  values  of  this  constraint  correspond  to  a  loss  of  information.  This  algorithm 
effectively  iterates  between  projections  in  the  object  plane  and  back  to  the  Fourier 
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plane  to  determine  the  CTHIS  reconstruction  similar  to  the  popular  Gercheberg- 
Saxton  algorithm  which  is  another  type  of  POCS  algorithm  [9]. 


1.6.3  Maximum  Likelihood  Estimators.  In  [7]  Descour  and  Dereniak  simply 
propose  mapping  each  volumetric  element  (voxel)  from  the  data  cube  to  the  pixel 
projection  in  the  image  plane.  They  experimentally  determined  the  system  transfer 
matrix  H: 

d(x,  y )  =  h(x,  y,  A)  *  o(x,  y,  A)  +  n(x,  y)  (1.18) 

where  d(x,y)  is  the  sensed  data,  o(x,y,  A)  is  the  spatial-spectral  object,  and  n(x,y) 
is  a  noise  vector.  Equation  (1.18)  directly  maps  the  discrete  image  pixels  to  their 
discrete  object  voxel  counterparts.  Descour  et  al.  [7]  chose  to  use  a  monochrometer 
and  physically  translated  the  input  at  each  (x,y)  value  determining  the  point  spread 
function  (. h(x,y ,  A))  directly  (each  voxel  to  pixel  projection  forms  a  single  PSF).  This 
step  maximizes  the  probability  density  of  the  object  given  the  data  Pr(o\d)  iteratively 
(equation  (1.19))  with  a  stopping  criterion  such  that  the  quotients  go  to  unity,  or  when 
|| d(x,  y)  —  h(x,  y,  A)  *  o(x,  y,  A)  1 1  is  minimized  below  a  certain  threshold. 


M 

ok+1(x,y,  A)  =  ok(x,y,  A)^  hm{x,y,  A)  * 

m=  1 


dm(x,y ) 

dkm(x,y) 


(1.19) 


In  this  step  d^x^y)  is  the  estimate  of  the  noisy  data  using  the  object  calculated  in 
the  previous  iteration. 


1.7  Reconstruction  Problems 

Even  with  algorithms  for  image  reconstruction  understood,  there  are  still  some 
problems  which  make  it  difficult  to  completely  restore  all  of  the  information  about 
the  object.  The  fact  that  the  sensor  cannot  fully  control  the  orientation  of  the  object 
leads  to  a  fixed  angle  of  projection.  The  result  is  that  the  system  transfer  matrix  has 
a  null-space  and  is  not  fully  invertible.  This  yields  two  related  problems,  first,  that 
we  desire  more  outputs  than  we  have  inputs  for  (the  result  is  that  the  problem  is 
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underdetermined),  and  second,  that  because  of  the  fixed  angle  of  projection,  we  have 
a  cone  of  missing  information  in  the  spatial-spectral  Fourier  domain  corresponding  to 
the  null-space  of  the  transfer  matrix. 

1.7.1  Underdetermined  Problem.  Given  either  CTH1S  design  discussed  pre¬ 
viously,  both  spatial  and  spectral  information  are  collected  simultaneously.  Because 
of  the  simultaneous  collection,  it  is  not  possible  to  directly  design  the  sensor  to  de¬ 
couple  the  spatial  and  spectral  resolutions.  Because  of  this,  changing  the  design  to 
increase  resolution  in  one  domain  (say  spatial)  affects  the  resolution  in  the  other  do¬ 
main.  Also,  because  of  the  fixed  angle  of  the  object  with  the  imager,  some  of  spatial 
frequencies  of  the  original  object  are  never  sampled.  The  transfer  matrix  for  this  case 
therefore  cannot  be  directly  inverted.  The  pseudo-inverse  of  this  matrix  (as  observed 
in  section  1.6.2),  will  provide  a  good  first  step  to  an  estimation  problem,  but  is  almost 
unusable  as  a  complete  object  estimate. 

1.7.2  Cone  of  Missing  Information.  As  mentioned  earlier,  a  CTH1S  image 
in  the  FPA  can  be  considered  as  a  projection  of  the  3D  data  cube  into  the  2D  CTH1S 
image  plane.  Figure  1.15  (b)  shows  the  data  cube  projected  onto  the  CTHIS  image 
plane.  The  angle  of  the  projection  is  determined  by  the  angular  deviation  from  the 
dispersion  element  for  the  wavelength  of  interest.  According  to  [2]  “the  2D  (Fourier) 
transform  of  a  2D  projection  yields  one  plane  through  the  3D  transform  of  the  original 
object”,  in  this  case,  the  original  object  is  the  hyperspectral  data  cube.  Each  of  these 
planes  pass  through  the  origin  and  the  angle  of  the  Fourier  plane  slice  is  orthogonal  to 
the  direction  of  the  corresponding  Fourier  axis  of  the  projection  from  the  data  cube 
onto  the  2D  CTHIS  image.  This  indicates  that  we  can  reconstruct  the  3D  Fourier 
data  cube  simply  by  taking  the  Fourier  transform  of  a  continuum  of  2D  projections 
(see  figure  1.17).  We  can  then  recover  the  hyperspectral  data  cube  simply  by  taking 
the  inverse  3D  Fourier  transform.  According  to  [7]  “Recovery  of  a  3D  distribution 
from  2D  projections  is  known  as  the  x-ray  transform.”  Also  according  to  [7],  the 
x-ray  transform  has  4  assumptions:  the  imager  is  continuous,  there  are  no  diffraction 
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effects  or  aberrations  (aberration  free  geometric  optics),  the  data  is  not  corrupted  by 
noise  (perfect  detector),  and  that  the  projections  can  be  obtained  at  any  azimuth  and 
projection  angles.  All  of  these  assumptions  are  violated  by  a  CTHIS,  especially  in  the 
presence  of  atmosphere. 

If  a  projection  of  the  data  cube  could  be  taken  at  any  azimuth  and  projection 
angles,  we  could  choose  a  projection  angle  perpendicular  to  the  wavelength  axis. 
Then  the  Fourier  plane  slices  are  perpendicular  to  the  spatial  frequency  axes  and 
then  we  can  simply  sweep  through  180  degrees  to  fully  reconstruct  the  3D  Fourier 
cube.  Unfortunately,  we  cannot  do  this  for  two  reasons,  first  the  projections  are 
taken  at  discrete  azimuth  angles  (finite  number)  which  yields  a  discrete  covering  set 
spanning  the  3D  Fourier  space,  and  second  we  cannot  take  a  90  degree  projection 
to  the  wavelength  axis.  The  finite  number  of  angles  can  be  overcome  by  accepting 
a  discrete  reconstruction  (which  we  have  already  accepted  from  any  hyperspectral 
data  cube),  in  other  words,  making  the  object  cube  discrete.  The  non-orthogonal 
projection  angle,  on  the  other  hand  will  yield  a  “cone  of  missing  information”  in  the 
Fourier  plane. 

This  cone  means  that  a  truly  unique  reconstruction  is  not  possible  from  the  mea¬ 
surements  alone,  however  a-priori  (object  or  transform)  domain  specific  constraints 
can  be  used  to  fill  in  the  missing  cone.  The  missing  cone  by  itself  shows  that  object 
features  which  have  higher  spatial-frequency  have  more  spectral  frequency  content 
represented  in  the  data  [7,20].  The  problem  of  reconstructing  the  cone  of  missing  in¬ 
formation  is  equivalent  to  the  limited-angle  computed-tomography  problem  [7].  Due 
to  the  cone  of  missing  information  and  the  fact  that  spatial  and  spectral  information 
is  correlated,  there  is  ambiguity  in  determining  the  resolution  of  these  sensors.  To 
date  no  study  of  CTHIS  resolution  has  been  performed. 

1.8  Literature  Study  Conclusions 

This  section  detailed  the  current  research  into  imaging  spectrometers  based  on 
computed  tomography.  Various  designs  of  other  hyperspectral  imaging  sensor  designs 
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were  discussed  as  well  as  some  of  the  advantages  and  disadvantages  of  each.  CTHIS 
was  discussed  starting  with  early  work  in  the  held  to  current  designs  using  diffractive 
elements,  prisms  and  chromatic  aberration  of  a  lens  for  providing  spectral  diversity. 
These  were  examined  with  advantages  and  disadvantages  from  each  one.  From  there, 
various  reconstruction  algorithms  were  discussed.  Finally,  the  problem  of  missing 
information  from  the  data  was  discussed  as  this  problem  affects  every  CTHIS  sensor. 
The  main  advantage  of  CTHIS  over  other  sensors  is  that  they  more  fully  utilize  the 
light  coming  in  from  a  scene  versus  diffractive  imaging  spectrometers  or  FTIS.  They 
also  have  the  ability  to  be  configured  to  image  hash  events  (unlike  a  diffractive  sensor), 
and  are  not  as  sensitive  to  vibration  due  to  their  simple  design  (unlike  a  FTIS). 
While  they  have  some  significant  advantages,  CTHIS  sensors  have  an  ambiguity  due 
to  the  general  nature  of  the  sensor  setups  with  the  cone  of  missing  information,  and 
this  leads  to  ambiguity  of  the  resolution  of  the  CTHIS  information.  Specifically  the 
spectral  resolution  has  not  been  studied  in  the  presence  of  noise.  Also,  very  little  has 
been  done  to  make  the  CTHIS  sensor  available  outside  the  laboratory  environment 
with  only  one  author  actually  referencing  data  taken  outside.  No  discussion  is  present 
in  the  literature  of  using  a  CTHIS  in  the  presence  of  atmospheric  turbulence  which 
would  be  necessary  for  the  use  of  CTHIS  in  real-world  applications.  Techniques  for 
evaluating  the  limiting  factors  of  the  CTHIS  spectral  resolution  and  solutions  for 
determining  the  spectral  resolution  will  be  discussed  in  the  following  chapters. 


II.  Theoretical  Lower  Bound  on  Resolution  of  CTHIS  in  the 

Presence  of  Noise 

This  chapter  develops  the  theoretical  lower  bound  on  the  spatial  and  spectral  reso¬ 
lution  of  an  unbiased  object  estimator  for  a  lens-based  CTHIS.  First,  the  variance 
obtained  by  the  CRLB  is  related  to  the  Rayleigh  resolution  criteria  in  section  2.1. 
Next,  a  description  of  the  Cramer-Rao  Inequality  is  developed  in  2.2.  Third,  an  im¬ 
age  model  is  developed  from  a  simple  object  model.  Finally,  the  lower  bound  on 
estimator  performance  for  a  simple  object  is  given  in  section  2.4. 

2. 1  Relationship  of  Estimator  Uncertainty  to  the  Rayleigh  Criteria 

The  Cramer  Rao  Lower  Bound  (also  called  the  CRLB)  gives  us  a  lower  bound 
on  the  variance  of  an  estimator.  However,  we  wish  to  determine  the  resolution  of 
our  system.  The  Airy  disc  is  the  PSF  from  a  circular  aperture.  LIsing  a  standard 
resolution  criteria  such  as  a  Rayleigh  criteria  [11],  two  point  sources  are  considered 
resolved  when  they  are  separated  so  that  the  peak  of  an  Airy  disc  from  one  point 
source  is  in  the  null  the  Airy  disc  of  another  point  source  next  to  it.  The  Airy  disc  is 
the  PSF  from  a  circular  aperture  which  has  a  peak  in  the  center,  and  concentric  rings 
of  reduced  intensity  with  null  rings  in  between.  In  the  case  of  a  circular  aperture, 
Figure  2.1  shows  an  example  of  two  points  barely  resolvable  with  the  Rayleigh  criteria. 
However  we  seek  to  apply  the  CRLB  in  an  attempt  to  capture  the  effects  of  noise  on 
resolution.  Specifically  when  the  actual  point  separation  is  equal  to  twice  the  standard 
deviation  of  the  estimate  of  the  separation,  we  can  say  that  we  are  statistically  resolved 
according  to  the  Rayleigh  criteria  of  separation  because  a  single  standard  deviation 
represents  the  uncertainty  of  points  in  either  spatial  or  spectral  dimensions.  So  the 
resolution  criteria  can  be  stated  as: 


A 

- =  1 

2o-a 


(2.1) 
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Figure  2.1:  Two  points  resolved  according  to  the  Rayleigh  criteria 

where  A  is  the  separation  of  the  two  points  in  question,  and  a  a  is  the  lower  bound 
on  the  standard  deviation  of  the  estimate  of  A  according  to  the  CRLB. 


2.2  Cramer- Rao  Lower  Bound  and  Fisher  Information  of  an  Image 
with  Poisson  Noise 

The  Cramer-Rao  inequality  states  that  the  variance  of  an  unbiased  estimator 
of  a  parameter  is  no  smaller  than  the  inverse  of  the  Fisher  Information  of  the  pa¬ 
rameter  estimated  [32] .  This  provides  a  lower-bound  on  the  variance  of  an  estimator, 
commonly  called  CRLB.  Specifically: 


>  FT1  =  -E 


d2L 

~d^ 


-i 


(2.2) 


where  a2  is  the  variance  of  the  estimated  parameter  0.  Also,  Lk  is  the  natural 
logarithm  of  the  conditional  Probability  Density  Function  (PDF)  of  the  data  with 
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respect  to  the  estimated  variable: 


Lk  =  \nP(dk\</>) 


(2.3) 


where  P(dk\(p)  is  the  probability  mass  function  of  the  data  d k  given  the  parameter 
(f).  If  more  than  one  parameter  is  being  estimated,  then  the  joint  probability  is  used, 
and  the  Fisher  Information  becomes  a  matrix.  For  multiple  parameters  estimation, 
we  refer  to  them  as  (pi  and  (p3  where  i  and  j  where  the  parameters  are  the  same  only 
when  i  =  j.  The  log-likelihood  over  all  data  collected  is  given  by  L.  Each  element 
of  the  matrix  is  calculated  from  the  second  partial  derivative  of  the  log-likelihood 
function  with  respect  to  the  parameters  being  estimated. 


y 


—E 


'  82L  ' 

_dotdOj_ 


(2.4) 


The  CRLB  for  each  parameter  is  the  main  diagonal  of  the  inverse  of  the  Fisher 
information  matrix  (i  =  j). 


>  F; 


-l 


—  *7 


(2.5) 


The  Probability  Mass  Function  (PMF)  of  a  collection  of  two  dimensional  images  with 
Poisson  noise  is  given  by: 


P(d\oi,  o2,  Ax,  AA)  =  nnn  ik(x,y)dk{x,y)e  lk{x’y)(dk(x,y)\ )  1  (2.6) 

y  x  k 

where  the  dk  indicates  a  single  frame  in  a  series  of  images.  The  intensities  of  two 
points  are  o\  and  02.  This  equation  assumes  that  the  location  of  the  first  point  is 
known  Only  two  points  are  used  because  this  will  yield  the  simplest  object  to  study 
both  spatial  and  spectral  resolution.  More  complex  objects  and  features  could  be 
estimated,  however  the  CRLB  requires  a  specific  object,  therefore  the  simplest  object 
is  chosen  for  this  resolution  study.  The  two  points  are  separated  by  (Ax,Aa).  The 
location  of  the  Erst  point  is  assumed  known  and  is  set  at  the  origin  for  convenience. 
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A  simple  object  is  given  by 


o(u,v,  A)  =  o15(A-A0)5('u-x0,'i;-7/o)+025(A-(Ao+AA))5('u^(xo+A;I.),n-7/o)  (2.7) 


With  this  simple  object,  the  resolution  of  imaging  systems  can  be  determined  by 
moving  one  point  with  respect  to  the  other,  and  using  the  resulting  image  to  compute 
the  CRLB.  The  image  equation  is  given  by: 

ik(x,  y)  =  y,  A)  +  o2hk(x  -  A x,y-  Ay,  A  +  Aa)  (2.8) 

A 

where  hk(x,y,  A)  is  the  point  spread  function  of  our  imaging  system.  As  mentioned 
in  section  2.1,  these  parameters  are  related  to  the  resolution  of  two  points  using  the 
standard  deviation  as  a  measure  of  the  ability  to  resolve  two  points.  This  relationship 
allows  a  determination  of  resolution  based  on  the  CRLB  and  gives  a  parameter  that 
can  be  used  to  predict  how  many  images  will  be  required  to  achieve  a  particular 
resolution.  The  log-likelihood  function  then  becomes: 


L  =  lnP(d\oi,  o2,  Ax,  AA)  =  EEE  dk(x,  y)ln[ik(x,  y)}  -  ik(x,  y)  -  ln[dk(x,  y)\] 

y  x  k 

(2.9) 

We  wish  to  find  the  lower  bound  for  the  variance  of  a  parameter  estimated  from  the 
data.  dk.  The  parameters  of  interest  will  be  discussed  in  the  next  section,  however 
we  can  greatly  simplify  the  CRLB  for  the  Poisson  PMF  with  generic  parameters. 
Computing  the  first  partial  derivative  with  respect  to  a  generic  parameter  of  interest: 


dL  y^  Y^  Y-^  (  7  /  \  1  dik(x,y)  dik(x,y)\ 

- 


(2.10) 


where  <f>i  is  the  parameter  of  interest.  Because  we  are  interested  in  more  than  a 
single  parameter,  the  second  partial  derivative  needs  to  be  taken.  The  second  partial 
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derivative  of  a  different  parameter  <f)j  is: 


d  L 

dcbidcbj 


= S  d*(x>  y) 

y  x  d 

dik(x:y )  dik(x,y)\ 


1  fdik(x,y) 
ik(x,y)2  \  d&dfij  %k  X,V 

d2ik(x,y) 


d<pi  d(j)j  J  dfpidcpj 


(2.11) 


Then  we  take  the  expected  value  to  get  the  Fisher  information  matrix,  keeping  in 
mind  that  the  expected  value  is: 


E[ik(x,y)]=  it(x,  y ) 


(2.12) 


that  is,  the  expected  value  of  the  noisy  data.  dk(x,  y)  is  the  image  ik(x,  y).  Using  this, 
the  second  derivative: 


-E 


d  L 


- E[dk(x,y )] 

y  x  d 

dik(x,y)  d ik(x,y)\ 


1  fd  2ik(x,y) 

ik(x,y)2\  d^d^j 

d2ik{x,y ) 


h{x,y) 


d(j>i  d(j)j  J  d  fad  fa 


(2.13) 


becomes  greatly  simplified: 


-E 


d  L 

dfadfa 


=  V  W  ( d2tk(x,y) 

ik(x,y)2{  dfadfa  k[  'V) 
_  d ik(x,y)  d ik(x,y)\  _  d2ik(x,y) 

dcj)i  dfa  J  dcf)id(j)j 


(2.14) 


And  after  cancelling  all  the  extra  terms,  the  final  form  of  the  Fisher  information 
matrix  elements  are  given  by: 


p  _  1  dik(x,y)  dik(x,y) 

13  yxk^^  d(t>i  ^ 


meaning  that  we  only  need  to  compute  the  first  derivatives  of  any  parameter  of  interest 
and  use  the  cross  products  to  compute  the  Fisher  information  matrix  elements. 
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2. 3  Image  model 


As  mentioned  above,  we  wish  to  find  a  bound  on  the  spatial  and  spectral  res¬ 
olution  of  our  CTHIS  system.  One  method  is  to  compute  the  estimator  variance 
for  a  scene  with  a  known  spectral  and  spatial  separation  (Aa,  A^),  and  compare  the 
standard  deviation  to  the  actual  separation.  The  scene  consists  of  two  point  sources 
separated  by  a  shift  of  (Aa,  Ax)  in  space  and  wavelength  respectively.  The  separation 
for  spatial  dimension  is  only  a  single  dimension  used  for  simplicity  (figure  2.2).  An 


Original  Object,  A  =5,A.=80nm 


Figure  2.2:  Example  Object 

image  of  a  simple  object  with  two  point  sources  at  (x0,  yo,  A0)  with  a  shift  of  (Ax,  Aa) 
between  them,  with  an  arbitrary  point  spread  function  (PSF)  is  modelled  by: 

h(x,y)  =££mo  -  A 0)6(u  -  x0,v  -  y0 ) 

A  u,v 

+  o2<5(A  -  (A0  +  AA ))5(u  -  (x0  +  Ax),v  -  r/0)]  ('2'16') 

x  hkiu  —  x.  v  — y.  A)  —  8 

where  k  is  an  index  to  a  particular  image  or  PSF  in  a  set.  The  variable  f3  is  a 
background  that  is  not  affected  by  the  PSF,  and  is  discussed  later.  Assuming  a 
spatially  shift- invariant  PSF  and  the  using  sifting  property  of  the  dirac-delta  function, 
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the  image  becomes: 


ik(x,y)  =  ^  Oi<S(A  -  A 0)hk(x  -x0,y-  y0,  A)- 

a  (2.17) 

+  o2S( A  -  (A0  -  Ax))  x  hk(x  -  (x0  +  A x),y  -  y0,  X)d,\  +  f3 

noting  that  the  Dirac  delta  for  the  spatial  dimension  of  the  point  source  fixes  the 
PSF  at  its  location.  Finally  defining  the  central  wavelength  Ao  and  integrating  over 
all  wavelengths  applies  the  sifting  property  of  the  Dirac  delta  yielding: 

ik(x,y)  =  oihk(x  -  x0,y  -  y0,  A0)  +  o2hk(x  -  (x0  +  A x),y  -  y0,  A0  +  Aa)  +  /3.  (2.18) 

Given  the  equations  developed  above  for  the  CRLB  (equation  (2.15)),  we  can  use  our 
image  model  developed  in  this  section  to  derive  a  more  specific  form  of  the  CRLB  for 
our  simple  test  object.  We  are  interested  in  our  ability  to  estimate  the  intensities  (oi 
and  o2)  and  the  relative  position  (Ax,  Aa)  of  the  two  point  sources  in  the  scene.  These 
are  the  specific  parameters  of  interest  used  for  fa  and  (p-j  in  the  previous  section  to 
compute  the  CRLB.  Notice  that  the  background  (3  is  a  constant  with  respect  to  the 
parameters  of  interest,  and  drops  out  after  taking  the  derivative.  It  will  only  affect 
the  divisor  from  equation  (2.15). 

dik(x,  y)  =  x  _Xo  y_  yo  ^q)  (2.19) 

GO  i 

'  i  '  ^  =  hk(x  —  (xq  +  A x),y  —  y0,  Ao  +  Aa)  (2.20) 

00-2 

Note  that  the  intensity  partial  derivatives  are  simply  the  fixed  point  spread  functions 
which  can  be  computed  directly.  However,  for  the  other  parameters  (Ax,  Aa) 

^  =  o2T^—hk(x  -  (x0  +  A x),y-  y0,  A0  +  Aa)  (2.21) 

a  Ax  a  Ax 

=  02J^hk(x  _  (Xo  +  A x),y  -  y0,  A0  +  Aa)  (2.22) 

(7Aa  GLA\ 
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more  simplification  will  be  needed. 


2.f  Theoretical  Lower  Bound  on  Spectral  Resolution  of  a  Defocused 
Image  with  Poisson  Noise 

We  need  to  determine  the  PSF  to  further  simplify  equations  (2.19)-(2.22)  in 
order  to  compute  the  CRLB.  The  intensity  point  spread  function  ( hk(x,y ,  A)  in  equa¬ 
tions  (2.21)  and  (2.22))  of  a  given  optical  system  can  be  computed  by  the  amplitude 
point  spread  function  h%(x,y,  A).  The  amplitude  point  spread  function  is  the  field 
resulting  from  a  plane  wave  propagating  through  the  optical  system,  corresponding 
to  a  point  source  (impulse)  at  infinity  (equivalent  to  the  impulse  response  discussed 
in  classical  linear  systems): 


hk(x,y,  A) 


\K(xix)\2 

L,yhl(x,y,  A) 


(2.23) 


Because  this  is  a  shift  invariant  amplitude  PSF,  the  denominator  will  later  be  shown 
to  be  simply  a  scaling  constant  which  does  not  vary  with  respect  to  x  or  A. 

In  the  case  of  a  defocused  lens  with  a  square  aperture  of  side  2w  and  a  defocus 
in  meters  of  Wd,  ignoring  scaling  constants  and  pure  phase  factors,  the  amplitude 
point  spread  function  (impulse  response)  is  computed  by  [11]: 


K(x,y,  A) 


rect{2wtf)rect{2wQ  exp 


-27rj  / Wd(£2  +  C2)  +  a[£H-j/C\ 

A  \  w2  zd  ) 


dfd( 


(2.24) 

where  (£,  ()  are  the  spatial  coordinates  corresponding  to  the  field  in  the  lens  plane, 
and  zd  is  the  distance  between  the  lens  and  the  focal  plane  for  frame  number  k.  Wd 
is  given  by: 


Wd  =  - 


(2.25) 


where  z%  is  the  distance  from  the  lens  to  the  image  forming  plane  (where  an  image 
would  be  in  focus).  The  variable  z%  depends  on  focal  length  which  is  dependent  on 
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the  index  of  refraction  (and  thereby  wavelength).  The  actual  focal  length  can  be 
calculated  from  material  properties,  however  for  the  sake  of  this  analysis,  the  change 
in  focal  length  due  to  a  change  in  wavelength  is  assumed  to  be  directly  proportional 
to  a  parameter  a,  with  a  possible  bias  term  (/3fi  a  minimal  focal  length)  therefore: 

/(A)  ccA  +  /3fi  (2.26) 

and  because  we  are  considering  a  point  source  at  infinity  (a  plane  wave  field  propa¬ 
gation  through  the  lens),  z,:  becomes: 

Zi{  A)  =  /(A)  (2.27) 


Therefore,  Wd  becomes: 


Wd 


w2  (  1  1  A 

2  \  za  aX  +  fifi J 


(2.28) 


Substituting  equation  (2.28)  into  (2.24)  we  get: 


K(x,y,X,d ) 


p  pw 

1  1  ovri 

~nj 

f(£2  +  C2) 

0 e2  +  c 2)  2(^+?/c)ai 

/  /  eXP 

J  J  —w 

X  1 

V  Zd 

aX  +  (3fi  zd  J  _ 

d^K 


(2.29) 


Where  the  limits  are  based  on  the  size  of  the  aperture.  Note  that  this  can  be  extended 
to  a  generic  aperture  function: 


K(x,y,X,d )  = 


^(f ,  0  exp 


nj((e  +  C2)  (^2  +  C2)  2(^  +  0 

A 


Zd 


aX  +  /3  f  i 


zd 


dt,d( 

(2.30) 


We  extend  this  to  a  generic  aperture  function  so  that  the  limits  of  our  equation  extend 
from  —  oo  to  oo  and  correspond  to  a  Fourier  transform. 


Because  we  are  also  estimating  the  intensities  of  the  point  source,  we  also  com¬ 
pute  the  Fisher  information  matrix  with  respect  to  the  parameters  Oi  and  02,  which 
are  measured  in  photons.  Because  the  intensities  in  number  of  photons  are  much 
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higher  quantities  than  the  pixel  size  and  the  wavelength  (typically  measured  in  mi¬ 
crometers),  the  Fisher  information  matrix  quickly  becomes  singular  as  the  number  of 
images  k  increases.  In  order  to  avoid  this,  we  make  the  parameters  a  discrete  step  size 
and  take  derivatives  with  respect  to  the  number  of  steps.  We  also  scale  our  aperture 
variable  to  keep  constant  sampling  in  the  focal  plane.  We  define  Ld  as  the  size  of  the 
aperture,  and  keep  the  ratio  fixed.  Keeping  this  ratio  fixed  scales  the  sampling 
to  a  constant  focal  plane  array  pixel  size  yielding  a  more  accurate  model  of  the  focal 
plane.  So: 

,  'Wei  , 

ax  —  T  ay  —  T 

Ld,  Ld 

.  1-j  d  *  -A  d 

At  =—  Ac  =— 

?  N  c  N  (2.31) 

x  —kdx  y  =pdx 

£  =™A?  C  =nA^ 

where  N  is  the  number  of  points  in  the  focal  plane.  We  assume  square  pixels,  and 
square  sampling  of  the  aperture  plane,  where  dx  is  the  pixel  size  in  the  focal  plane, 
and  A^  is  the  sampling  in  the  aperture  plane.  Now  to  compute  the  normalized  Point 
Spread  Function,  we  substitute  equations  (2.30)  and  (2.31)  into  (2.23),  and  replace 
the  integrals  with  summations. 


hk(k,p,  A)  = 


A|  E  A(m,  n)  exp 

m,n  L 

Af  ^2  A(m,n)ex p  — 

m  n 


t rj(m2+n2)  f  A 2zd _ A-i-d 

A  \^N2dx2  ( a\+(3fi)N2dx 2 


■rrj(in2+n2)  (  A 2zd _ A2"d 

A  \N2dx 2  (aX  -\-/3fi)N2dx2 


j  exp  — 
E  exp 

T  'll 


—2nj  (km-\-np) 
N 


—27 vj  ( km-\-np ) 
N 

(2.32) 


As  seen  in  equation  (2.15),  the  Fisher  information  matrix  elements  are  computed  by 
the  first  derivatives  of  the  image  with  respect  to  the  variables  under  consideration. 


Specifically  we  are  interested  in  the  Ax  and  AA  paramaters.  Using  a  discrete  change 
of  variable  for  Ax  and  A^.  The  variables  A'x  and  A^  signify  a  step  size  in  meters  (for 
simplicity,  we  fix  A^.  as  the  size  of  1  pixel),  and  Px  and  P\  signify  the  number  of  steps 
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in  discrete  space: 


Ax  =  A  'XPX=  dxPx 
Aa=  A  xPx 

So  the  the  image  equation  (2.34)  becomes: 


(2.33) 


h(k,p)  =  o1hk(k,p,  \0)  +  o2hk(k- Px,p,\0  + A'XPX)  (2.34) 


assuming  xq  =  0  and  yo  —  0  for  the  sake  of  notational  simplicity.  The  derivatives  are 
now  taken  with  respect  to  Px  and  P\.  This  necessitates  taking  the  derviatives  of  the 
PSF  with  respect  to  these  variables.  Utilizing  the  product  rule  and  an  identity  of  the 
absolute  value  function,  the  derivative  of  the  PSF  with  respect  to  Px  becomes: 


d  d 

hk(k  -  Px,p,  A0  +  A 'XPX)  =— -K~lhak{k  -  Px,p,  A0  +  A XPX) 


OP; 


3PX 

x  hf{k-px,p,  A0  + A^Pa) 
d 


=2k  1Real 


dPr 


hak(k-Px,p,\0  +  AxPx) 


(2.35) 


x  hk*(k  —  Px,p,  A0  +  A'xPx) 


leading  to: 


=2  k  1Real 


ST  a  /  \  2njm 

2_^  A{m,  n)  —  exp 


N 


nj(m2  +  n2)(A0  +  A  xPx)zj 
N2dx 2 


x - 


Zd  (a(Ao  +  AaPa)  +  j3fi)  J  \ 


exp 


—27 \j{km  +  np) 


exp 


2njrnPx 

N 


x 


exp 


exp 


N 

A(m,  n) 

m,n 

nj{m2  +  n2)(A0  +  A'xPx)z^  (  1 


N2dx2 

—2it j[km  +  np) 
N 


exp 


Zd  (a(Ao  +  A  xPx)+/3fi) 

2irjmPa 
TV 


(2.36) 
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where  k  is  a  scaling  constant  normalizing  the  area  under  the  PSF  to  unity  and  is  the 
same  as  the  denominator  in  (2.23).  The  derivative  of  the  PSF  with  respect  to  PA  is 
given  by: 


d 

dFx 


K(k 


PxiPi  A o  +  AaPa) 


1hl(k-Px,p,Xo  +  A'xPx) 

x  hl*(k  —  Px,p,  A0  +  AaPa) 


=2k  lReal 


d 

ap 


hak(k-Px,p,\0  +  A'xPx) 


(2.37) 


x  h%*(k  —  Px,p,  A0  +  AaPa) 


yielding: 


=2k  1Real 


A(<m' 


n) 


irj(m 2  +  n2)A'xz^ 


N2dx2 


Zd  (a(Ao  +  A  XP\)  +  (3fi) 
nj{m2  +  n2)(A0  +  A  'xP\)z2d 


«A, 


N2dx2 


(a(A0  +  A  aPa)  +  fdfi )2 


exp 


nj(m2  +  n2)(A0  +  A'xP\)zj 


1 

x  I  — 


TVW 

1 


Zd  (a(Ao  +  A  XP\)  +  fdfi) 


x  exp 


-27 Tj(km  +  np) 


TV 


exp 


2/njmPx 

TV 


x 


A(m,  n)  exp 


7rj(m2  +  n2)(A0  +  AaPa)^ 
N2dx2 


Zd  (a(Ao  +  A  aPa)  +  /3/z) 


exp 


-2nj{km  +  np) 
TV 


exp 


2njmPx 

TV 


(2.38) 


It  should  be  noted  that  the  discrete  summations  in  equations  (2.32),  (2.38)  and  (2.36) 
can  all  be  computed  by  using  a  Discrete  Fourier  Transform  (DFT),  and  as  such,  the 
DFT  is  conveniently  used  to  compute  these  functions. 
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III.  Data  simulation  and  reconstruction 


This  chapter  discusses  the  development  of  a  projection-based  Expectation  Maximiza¬ 
tion  (EM)  Estimator  based  on  Poisson  noise  in  section  3.1,  followed  by  a  discussion 
of  the  simulation  setup  in  section  3.2. 

3.1  Projection  Based  Reconstructor 

In  this  section,  a  reconstructor  to  estimate  the  object  from  a  series  of  images 
is  developed.  This  reconstructor  is  the  basis  for  the  simulation  and  the  laboratory 
experiment  results  shown  in  section  5.4.2.  The  lens-based  chromotomagraphic  hyper- 
spectral  sensor  takes  a  series  of  images  each  having  a  known  defocus  from  the  next 
image  in  the  series.  Given  the  collected  image  data,  the  originating  scene  is  unknown 
so  an  estimator  the  must  be  used.  A  well  known  technique,  known  as  maximum- 
likelihood  estimation  (MLE)  is  to  maximize  the  conditional  probability  over  all  values 
of  a  parameter  to  be  estimated.  In  this  case,  the  unknown  parameter  is  the  scene  and 
the  data  is  the  image  produced  by  our  known  transfer  function  [32], 

The  reconstructor  is  further  simplified  to  a  projection  based  reconstructor  as 
we  only  intend  to  explore  the  spatial  separation  in  one  dimension,  as  well  as  spectral 
separation.  The  data  is  modeled  as: 

dk(x,y)  =  ik{x,y)  +nk(x,y)  (3.1) 

and  is  a  random  variable  with  a  Poisson  probability  mass  function  given  by  Po(d ): 

pD (d)  m 

y  x  k 

where  ik(x,y)  is  the  expected  value  of  the  data  dk(x,y)  for  frame  number  k  with  the 
2  dimensional  location  (x,y),  this  assumes  that  each  pixel  is  statistically  independent 
from  the  others.  The  value  ik(x,y)  is  also  the  image  formed  by  the  unknown  object 
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at  plane  k.  Specifically,  the  image  is  given  by: 


ik(x,  =  7a  o(u,  v,  A  )hk(u  -  x,  v  -  y,  A)  (3.3) 

A  u  v 

where  7^  is  the  spectral  weighting  paramater  for  the  number  of  photons,  o(u,v,  A) 
is  object  at  spatial  coordinates  (u,  v)  in  the  object  plane  and  wavelength  A,  and 
hk(x,y,  A)  is  the  point  spread  function  at  frame  number  k.  The  spectral  weighting 
parameter  is  given  so  that  initially  °(ui  v->  x)  =  1  that  is  the  object  is  unity 

for  each  wavelength.  The  MLE  is  developed  by  maximizing  the  conditional  probability 
given  an  image  ik(x,y )  using  an  object  estimate  o(x,y,  A)  and  a  spectral  weighting 
estimate. 

Because  the  natural  logarithm  is  an  increasing  function,  the  computation  of  the 
probability  mass  function  can  be  greatly  simplified  by  applying  the  natural  logarithm 
and  maximizing  it  instead: 

LD(d)  =  ln[PD(d)]  =  ^2^2^2dk(x,y)ln[ik(x,y)\  -  ik{x,y )  -  ln[dk(x,y)\]  (3.4) 

y  x  k 

the  logarithm  of  the  factorial  term  in  this  equation  does  not  contribute  to  the  maxi¬ 
mization  and  is  ignored  in  further  development.  This  likelihood  function  is  recast  in 
terms  of  the  estimated  parameters. 

A  different  method  of  computing  the  object  was  proposed  in  [12]  where  a  pro¬ 
jection  operation  produced  a  single  spatial  dimension  and  a  spectral  dimension.  This 
method  features  a  significantly  lower  computation  cost.  The  image  is  related  to  a 
one-dimensional  image  projection  i^{x)  by  the  equation: 

h(x)  =  X)h*(x  ~u,y- v, x)  (3-5) 

y  A  u  v 
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Figure  3.1:  Description  of  the  EM  algorithm 
by  using  the  projection  of  the  object  and  PSF: 


ik(x)  =  EE  7ao(m,  X)hk(x  -  u,  A)  (3.6) 

A  u 

where  hk(x,  A)  =  hk(x,y,  A)  and  d(w,  A)  =  ^2vo(u,v,  A).  Using  these  projections 
significantly  reduces  computation  time,  but  still  allows  the  study  of  the  effect  of  the 
number  of  defocus  frames  on  spatial  and  spectral  resolution. 

The  next  step  in  estimating  the  scene  is  to  maximize  the  log-likelihood  which 
we  call  Q(o,  7a)  restated  in  terms  of  our  estimated  variables  with  respect  to  an  object 
estimate  o(u,  A)  and  spectral  weighting  estimate  7a-  This  could  be  done  in  multiple 
ways,  however  a  well  known  method  is  the  Expectation  Maximization  (EM)  algorithm 
[6].  In  order  to  use  the  EM  algorithm,  we  first  surmise  the  existence  of  an  a  set  of 
unknown  variables.  We  also  surmise  that  there  exists  a  mapping,  in  this  case  a  many- 
to-one  mapping,  between  the  set  of  unknown  variables  and  the  measured  data.  The 
EM  algorithm  then  maximizes  the  PMF  of  the  unknown  variables  conditioned  on 
the  measured  data  using  the  conditional  expected  value  derived  from  the  assumed 
mapping.  Figure  3.1  gives  a  flowchart  for  the  EM  algorithm.  Stated  simply,  we  seek 
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an  object  d(u,  A)  and  spectral  energy  parameter  7a  with  Q(dold,rfx°ld )  such  that: 


g(dne",  7Ane“’)  >  Q(dold,  rf\old)  (3.7) 

where  dnew  and  'yxnew  are  iterative  updates  to  our  previous  (old)  object  and  spectral 
energy  estimates  respectively.  We  begin  by  modelling  the  measured  (incomplete) 
data  dk(x )  as  a  sum  of  the  unknown  random  variables  dk(x)  =  Y^u  dk(x\u),  where  the 
collection  of  unknown  variables  dk{x\u)  is  called  the  complete  data.  In  creating  the 
complete  data,  each  variable  dk{x\u)  is  assumed  to  be  independent  for  each  unique 
value  of  u  and  Poisson  distributed.  The  complete  data  only  needs  to  be  statistically 
consistent  with  our  measurements  [27]  and  may  not  have  a  physical  meaning  as  is  the 
case  in  our  problem  statement.  The  mean  of  the  complete  data  is  given  by: 

E[dk{x\u)\  =  7a o(u,  A )hk(x  -  u,  A)  (3.8) 

The  complete  data  is  then  further  related  to  the  modelled  data  and  image  by: 

E[dk(x)]  =  ^2  E[dk(x\u)]  =  22  7a«(m,  A )hk(x  -  u,  A)  =  ik(x)  (3.9) 

U  U 

yielding  an  expression  from  which  we  can  develop  the  complete  data  log-likelihood. 
The  complete  data  log-likelihood  is  given  by: 

L(°,  7a)  =  ZEE  dk(x\u)ln[y\o(u,  A )hk(x—u,  A)]  —  EEE  7ao(w,  X)hk(x—u,  A) 

k  x  u,  k  x  u 

(3.10) 
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The  incomplete  data  log-likelihood  is  given  by  the  expected  value  of  the  incom¬ 
plete  log-likelihood  conditioned  on  the  measured  data  dj~{x)\ 

Q(6new,ixnew)  =E  [. L(6new,ixnew)\dk(x )] 

=  E  E  E  Eoiyk(x\u)\dk(x)]ln[ixnewdnew(u,  X )hk(x  -  u,  A)] 

k  x  u 

-  V  V  V  7-A~«8™"’(a,  A )hk(x  -  u,  A) 

k  u  x 

(3.11) 

where  the  expected  value  is  based  on  old  estimates  (designated  by  the  operation 
Eold[ *])  for  the  object  ( dold )  and  spectral  energy  coefficient  (rf\old).  We  designate 
the  object  (yonew)  and  spectral  energy  coefficient  ( yAne“’ )  as  iterative  updates  whose 
equations  are  derived  later,  the  old  values  of  which  are  previous  estimates  to  the 
data.  This  portion  of  the  process  where  we  develop  the  conditional  log-likelihood 
whose  expected  value  is  the  log-likelihood  of  interest  is  called  the  expectation  step 
of  the  EM  algorithm.  At  this  point,  it  is  necessary  to  develop  the  expression  for  the 
expected  value  of  a  single  Poisson  variable  conditioned  on  the  sum  of  a  number  of 
independent  Poisson  variables,  in  this  case,  dk(x\u)  conditioned  on  dk(x). 

E[dk(x\u)\dk(x)\  (3.12) 

In  order  to  determine  the  expected  value  of  a  single  Poisson  variable  conditioned  on 
the  sum,  we  start  with  the  sum  of  2  independent  Poisson  variables: 

d  =  di  +  d,2  (3.13) 


with  the  expected  values: 

E[di\  =  ii 
E[d2]  =  i2 


(3.14) 
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where  d\  is  the  Poisson  random  variable  of  interest  given  the  conditional  expectation 
of  the  poisson  incomplete  data  dk(x).  We  wish  to  find  the  PMF  of  d\  conditioned  on 
the  sum  d.  Using  Bayes  rule  we  get  the  following, 


P(Di\D)(di\d)  — 


P{Di,D2)){d\,d2) 

W) 


(3.15) 


that  is,  the  joint  density  of  d\  and  d2  divided  by  the  marginal  density  for  the  sum. 
Because  d\  and  d2  are  independent,  the  joint  density  is  given  the  product  of  their 
marginal  Poisson  mass  functions: 


P(Dl,D2)(di,d2)  —  Pdj,  {d'i)  Pd2  (^2)  — 


ifi22e"(il+*2) 


d\\d2\ 


(3.16) 


and  the  marginal  is  given  by: 


Pn(d)  = 


•  r]  — j 

1  e 


d\ 


(3.17) 


Noting  that  the  expected  value  of  d  is  the  sum  of  the  expected  values  of  d\  and  d2 


i  =  E[d ]  =  E[di\  +  E[d2\  —  i\  +  i2 


(3.18) 


and  also,  that  d2  can  be  re-written  as  a  difference  between  the  sum  and  d\ 


d2  =  d  —  d\ 


(3.19) 


we  can  find  the  conditional  expected  value.  Substituting  (3.19)  into  (3.16)  and  (3.18) 
into  equation  (3.17),  we  can  see  that  (3.15)  becomes: 


P{D1\D){di\d)  — 


d\ 


jdijid—d,!) 


1  °2 


d\ 


jdij{d—d\) 


1  "2 


di\(d  -  di)\  (ii  +  i2)d  di\(d  -  di))!  (U  +  i2)dl(U  +  u)d(U  +  u)_dl 

(3.20) 


we  can  rearrange  (3.20)  into  the  following  form. 


P(D1\D)(di\d)  — 


d\ 


i1  \dl  (  i2  ^  {d~dl) 


di\(d  -  di))!  \(*i  +i2)J  \(ii  +  i2) 


(3.21) 
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Noting  that  the  two  fractions  sum  to  unity: 


k  _  k  + k  k  _  k 
k+k  k  +  k  k  +  k  k  +  k 


(3.22) 


it  becomes  apparent  that  this  conditional  probability  is  in  the  form  of  a  binomial 
distribution.  A  binomial  distribution  in  general  form  is  given  by: 


PK(k) 


n\ 


k\(n  —  k)\ 


pk{l-pfn~k) 


(3.23) 


where  k  is  the  random  variable  out  of  n  independent  experiments  (called  Bernoulli 
trials)  with  a  probability  of  success  p.  The  expected  value  of  a  binomial  random 
variable  is  given  by  rip.  Letting  k  =  d\  as  our  random  variable,  with  n  —  d  trials  and 
also  letting: 


k 


P  =  — — 
ll  +  12 


(3.24) 


the  expected  value  of  a  single  Poisson  conditioned  on  the  sum  becomes  the  following 


£[di|d]  =  d- 


L  +  12 


(3.25) 


Substituting  in  dk(x)  and  k  =  E[dk(y\x)]  into  3.25  we  get  the  following(  [28]): 


E[dk(x\u)\dk(x)]  =  dk{x 


7ao(m,  X)hk(x-u,X) 


ik{x 


(3.26) 


It  is  also  desirable  to  incorporate  the  effect  of  background  noise  in  order  to 
make  the  simulation  more  realistic.  A  flat  Poisson  background  bk(x)  with  a  constant 
mean  (E[bk(x)]  =  c )  is  incorporated  by  adding  it  to  the  previous  data  model  dk(x,y) 
(equation  (3.9))  and  the  incomplete  data  is  relabeled  d\(x)  to  signify  the  addition  of  a 
background  dbk(x,  y )  =  dk(x,  y)  +  bk(x,  y).  This  background  is  similar  to  adding  a  flat 
held  to  the  image.  The  background  in  this  case  is  assumed  to  be  known  or  measured 
separately  and  is  not  estimated  in  this  algorithm.  This  model  is  effectively  diffuse 
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background  light  and  can  also  be  considered  similar  to  the  effect  of  additive  Gaussian 
noise  if  the  diffuse  light  level  is  high  enough. 

Since  we  are  only  interested  in  the  portion  of  the  data  dependant  on  the  object, 
bk{x)  can  be  considered  another  unknown  variable  added  together  with  the  complete 
data  with  a  total  expectation  of: 

E[dbk(x)}  =  ik(x)  +  bk(c )  (3.27) 


However,  because  we  measure  the  data  in  our  experiment,  the  measured  or  known 
background  is  used  instead  of  estimating  the  expected  value  and  no  further  estimator 
for  the  background  is  developed.  This  means  that  the  expected  value  in  (3.26)  is 
given  by: 


E°  ci[dk(x\u)\dbk(x)\  =  dk(x 


7 fddold(u,  X)hk(x  —  u,  A) 


(3.28) 


ik{x)  +  bk{x) 

where  ik{x)  is  the  image  formed  at  distance  plane  k  by  a  current  object  estimate  6old 
and  bk{x)  is  the  background  measured  previously. 


The  next  step  is  to  maximize  the  complete  data  log-likelihood  which  is  called 
the  maximization  step  of  the  EM  algorithm.  To  maximize  the  complete  data  log- 
likelihood,  the  derivatives  of  equation  (3.11)  are  taken  with  respect  to  dnew  and  'yxnew 
giving  equations  (3.29)  and  (3.30). 


dQ  _  Eold[dk(x\u)\dk(x) 

ddnew(u.X)  ~  2-^  2-^  dnew(u.  A) 

k  x 

Where  maxima  exist,  the  first  partial  derivative  will  equal  zero. 

dQ  _  Eold[dk(x\u)\dk(x)] 

Qrynew  /  V  /  V  /  v  rytiew 

^  k  x  u  ^ 

k  x  u 


(3.29) 


(3.30) 
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We  then  set  these  equations  equal  to  zero  solve  for  dnew  and  7™ ew .  We  further  make 
the  assumptions: 


^  hk(x,  A)  =  1 


X 


(3.31) 


it 


which  can  be  controlled  based  on  our  initial  guess  for  the  object  and  by  scaling  the 
PSFs.  The  spectral  energy  parameter  contains  the  energy  in  the  series  of  data  frames, 
and  is  assumed  equal  over  all  wavelengths  to  start  with: 


^ old  _  12  k  12  x  d'k(x) 


(3.32) 


where  e  is  the  number  of  wavelengths  of  interest.  Therefore  the  spectral  energy  is 
weighted  equally  across  all  wavelengths  according  to  the  total  energy  in  the  collection 
of  frames.  Also,  we  assume  that  for  each  iteration 


7^e“’  ps 


if 


(3.33) 


meaning  that  we  do  not  expect  the  estimate  for  7^  to  vary  significantly  with  each 
iteration.  Finally,  we  need  to  take  the  expected  value  of  the  incomplete  data  condi¬ 
tioned  on  the  complete  data  with  background  given  in  equation  (3.28).  Substituting 
the  expected  value  (3.26)  into  equations  (3.29)  and  (3.30)  gives: 


k  x 


d\{x) r)\doold(u,  A )hk{x  —  u,  A) 
(ik(x)  +  bk(x))dnew(u,  A) 


-7»“”  A) 

k  X 


0 


(3.34) 


k  X  u 


dbk(x) ryfdold(u,  A )hk(x  —  u,  A) 
(ik(x)  +  bk(x))^lew 


k  u  x 

(3.35) 
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with  the  summations  rearranged.  The  PSFs  are  shift  invariant  and  because  of  our 
earlier  assumption  of  Yhx  hk(x  —  u,  A)  =  1,  we  get: 


d\(x)ryfddold(u ,  A )hk(x  —  u,  A) 
V  *  ( (x)  +  h(x))dnew(u,  A) 


(3.36) 


SEE 

k  x  u 


d\(x)^fddold(u,  A )hk(x  —  u,  A) 
(ik(x)  +  6fc(a;))7))ew 


—  K 


0 


(3.37) 


where  the  value  K  is  the  total  number  of  frames.  Solving  the  equations  for  our 
estimates,  we  get: 


old 

dnew(u,  A)  =A'-14,o»mKA) 
7a 


EE 

k  x 


dbk(x ) 


ik(x )  +  bk{x) 


hk(x,  A) 


(3.38) 


7a”  =J 


r‘EE  jttwtIt)  E  A)M*,  A) 


(3.39) 


Finally  we  use  the  assumption  in  (3.33)  to  cancel  the  dependence  on  from  equation 
(3.38)  giving  us  the  final  form. 


40*0 


hfe(x,  A) 


(3.40) 


Equation  (3.40)  is  very  similar  in  form  to  the  object  expressed  in  [12],  and  the  added 
spectral  energy  term  (3.39)  allows  us  to  constrain  total  energy  in  the  object  estimate 
while  allowing  the  energy  in  each  wavelength  to  vary  independently.  This  solution  is 
similar  in  form  to  the  equations  in  section  1.6,  it  differs  on  two  key  points.  First,  7a 
is  allowed  to  range  and  is  estimated  jointly  with  the  object.  This  helps  to  separate 
the  spatial  and  the  spectral  dimensions  of  the  estimate.  Secondly,  this  algorithm  also 
has  an  independent  background  which  could  allow  this  algorithm  to  perform  better 
in  the  case  of  a  large  amount  of  diffuse  clutter. 
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3. 2  Simulation  Setup 

A  simulation  was  performed  using  the  parameters  listed  in  table  3.1.  Varying 
numbers  of  frames  and  spectral  separations  were  used  to  accurately  quantify  the 
relationship  between  spectral  resolution  and  the  number  of  frame  positions  (also  called 
defocus  planes). 

A  simple  experimental  setup  was  modelled  in  simulation  and  set  up  in  the 
laboratory  for  a  single  instance  of  the  simulation.  Figure  3.2  gives  a  notional  schematic 
for  the  optical  design  of  a  lens-based  CTH1S.  The  parameters  for  the  simulation  and 
experiment  setup  are  given  in  table  3.1.  Light-emitting  diodes  were  simulated  to 
approximate  the  object  modelled  in  the  CRLB  mentioned  in  section  2.4. 

The  source  object  was  modelled  with  a  circular  pattern  at  multiple  wavelength 
spacings  (Aa),  and  then  convolved  with  the  PSF  (over  multple  wavelengths)  corre¬ 
sponding  to  a  lens  with  chromatic  aberration. 

The  dispersion  element  and  the  two  lenses  closest  to  the  camera  were  simu¬ 
lated  as  a  single  converging  lens  for  this  setup,  and  the  dispersion  is  provided  by  the 
chromatic  aberration  of  the  converging  lens.  This  chromatic  aberration  changes  the 
focal  length  at  which  each  wavelength  gives  the  sharpest  picture.  The  wavelength  is 
modelled  as  in  the  CRLB  section  as  a  linear  change  given  by: 

fx  =  aX  +  fbias  (3-41) 


where  a  is  the  dispersion  parameter  measured  from  the  lens  and  A  is  the  wavelength. 
The  fbias  parameter  is  inherent  to  the  lens  and  is  dependent  on  both  the  radius  of 
curvature  and  the  index  of  refraction.  Both  a  and  fbias  were  measured  in  the  lab,  and 
used  in  simulation.  Equation  (3.42)  gives  the  un-normalized  PSF  of  the  lens  which  is 
normalized  to  unit  magnitude  for  the  reconstruction  algorithm  (see  previous  section). 


hk(x,y,  A) 


7T  j  f  ( U 2  +  V  2) 

A  V  zd 


(■ u 2  +  v 2) 

oA  +  fbias 


2  [xu  +  yv) 
Zd 


2 


dudv 


(3.42) 
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Figure  3.2:  Schematic  of  experiment  set  up  with  two  sources 

The  distance  Zd  represents  the  physical  distance  between  the  lens  and  the  detector 
plane.  Multiple  distances  for  Zd  are  used  by  moving  the  lens  further  away  from  to  the 
detector  plane  from  an  inital  starting  point.  The  initial  starting  point  is  determined 
by  the  wavelengths  of  interest  and  the  parameters  of  the  lens.  Specifically,  the  shortest 
distance  to  the  focal  plane  should  be  closer  than  the  focal  length  f\  of  the  shortest 
wavelength  of  interest.  Similarly,  the  farthest  Zd  from  the  focal  plane  should  be  farther 
away  than  the  focal  length  of  the  longest  wavelength  of  interest. 

The  number  of  defocus  frames  was  varied  giving  various  reconstructions  of  the 
object  to  determine  the  specific  effect  of  on  spectral  resolution.  A  multiplier  was  used 
to  vary  the  exposure  time  so  that  higher  numbers  of  images  did  not  result  in  larger 
numbers  of  photons  thereby  resulting  in  an  automatically  higher  signal  to  noise  ratio. 
The  exposure  time  was  varied  corresponding  to  the  number  of  frames  used,  with  a 
longer  exposure  (more  images  added  together)  corresponding  to  fewer  frames,  and 
shorter  exposure  (fewer  images  added  together)  corresponding  to  more  frames  used. 
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Table  3.1:  Simulation  and  Experiment  Parameters 
Parameter  Value  Units 


dx  (Detector  Pitch) 

16 

/irn 

N  (Detector  size) 

512 

pixels 

Aperture  Diameter 

1 

cm 

Number  of  wavelengths 

10 

Ao  (center  wavelength) 

560 

nm 

AA  (wavelength  spacing  simulation) 

10-100 

nm 

Minimum-Maximum  simulated  wavelengths  ( Amm  —  Amax) 

560-645 

nm 

fx0  (Focal  length  @  A0) 

441.6 

mm 

a  measured  from  lens) 

1.2406xl05 

fbias  (minimum  focal  length) 

373.4 

mm 

K  (Number  of  defocus  planes) 

3-20 

Minimum/Maximum  Distance  (zd) 

420.4/469.2 

mm 

AZd  (Distance  between  each  defocus  plane) 

0.26 

mm 

Number  of  iterations  for  reconstruction 

2000 

Number  of  noisy  realizations 

500 

This  normalizes  the  energy  (number  of  photons)  used  for  each  set  of  defocus  frames 
so  that  taking  more  frames  will  not  be  better  simply  because  of  a  longer  exposure 
time.  A  constant  background  was  then  simulated  and  then  a  noisy  realization  was 
taken  (using  Poisson  statistics).  The  simulated  images  and  background  were  scaled 
according  to  the  parameters  of  a  realistic  camera  to  convert  the  units  to  photons 
based  on  the  photon  transfer  characteristics  of  the  camera. 
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IV.  Experimental  Setup 

This  chapter  discusses  the  experimental  parameters  in  Section  4.1  and  the  reasons 
for  design  choices  for  the  experiment.  Then  in  Section  4.2  the  calibration  of  sources 
using  a  double-slit  experiment  is  discussed.  Finally,  modelling  calibration  of  the  lens 
and  the  comparison  to  the  CRLB  are  discussed  in  Section  4.3. 

4-1  Experimental  Parameters 

The  two  sources  were  two  light-emitting  diodes  (LED)  with  a  pinhole  in  front 
of  them  to  approximate  the  object  modelled  by  the  CRLB  in  the  previous  section. 
The  sources  were  aligned  using  a  beamsplitter  (to  combine  the  two  LEDs)  so  that 
they  overlapped  to  give  the  appearance  of  a  single  source.  After  the  beamsplitter, 
the  source  was  collimated  using  an  achromatic  collimating  lens,  so  that  the  distance 
between  the  source  and  dispersion  element  did  not  affect  the  measurement.  In  the 
experimental  setup,  only  two  diodes  with  center  wavelengths  at  at  560nm  and  645 nm 
were  used  in  the  setup.  The  center  wavelengths  were  determined  by  a  Young’s  double- 
slit  experiment.  The  double-slit  experiment  is  discussed  in  more  detail  in  section  4.2. 

The  fbias  parameter  is  inherent  to  the  lens  and  is  dependent  on  both  the  radius  of 
curvature  and  the  index  of  refraction.  Both  a  and  f^as  were  measured  in  the  lab,  and 
used  in  simulation.  Equation  (3.42)  gives  the  un-normalized  PSF  of  the  lens  which  is 
normalized  to  unit  magnitude  for  the  reconstruction  algorithm  (see  previous  section). 
Experimentally,  multiple  images  were  taken  at  each  defocus  plane.  The  images  were 
then  added  together  to  simulate  varying  exposure  times.  In  the  experimental  setup, 
there  was  a  significant  amount  of  stray  light  present  adding  to  a  non-wavelength 
specific  diffuse  background.  To  compensate  for  the  diffuse  lighting  experimentally, 
a  background  image  was  taken  along  with  each  defocus  plane.  If  a  background  is 
otherwise  unavailable,  the  amount  of  background  could  also  be  estimated  after  the 
fact  similar  to  the  derivation  given  above  for  estimating  the  image  from  an  object 
estimate.  However  this  research  simply  uses  a  measured  background. 
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Table  4.1:  Experiment  Parameters 

Parameter  Value  Units 


dx  (Detector  Pitch) 

16 

/jrn 

N  (Detector  size) 

512 

pixels 

Aperture  Diameter 

1 

cm 

Number  of  wavelengths 

10 

A0  (center  wavelength) 

560 

nm 

AA  (wavelength  spacing  simulation) 

10-100 

nm 

Minimum- Maximum  simulated  wavelengths  (A min  ~  A max) 

560-645 

nm 

Experimental  wavelengths 

560/645 

nm 

f\0  (Focal  length  @  A0) 

441.6 

mm 

a  (^  measured  from  lens) 

1.2406x10s 

fbias  (minimum  focal  length) 

373.4 

mm 

K  (Number  of  defocus  planes) 

3-20 

Minimum/ Maximum  Distance  (z^) 

420.4/469.2 

mm 

AZd  (Distance  between  each  defocus  plane) 

0.26 

mm 

Number  of  images  per  defocus  plane 

500 

Number  of  iterations  for  reconstruction 

2000 

Number  of  noisy  realizations 

500 

Diode  1  intensity  (550nm) 

3.64  x  1018 

Photons/ frame 

Diode  2  intensity  (645nm) 

1.82  x  1018 

Photons/frame 

Background  intensity  (645nm) 

3.56  x  1012 

Photons/frame 
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A  digital  camera  only  measures  light  intensity  in  digital  counts.  Because  the 
projection-based  reconstructor  assumes  Poisson  statistics,  the  measured  data  needs  to 
be  scaled  to  reflect  the  number  of  photons  sensed,  not  the  number  of  photoelectrons 
generated  (whole  number  digital  counts).  Knowing  that  the  Poisson  statistics  have 
both  an  expected  value  that  equals  the  variance,  the  data  was  assumed  scaled  so  that 
the  variance  matched  the  mean  by  dividing  by  the  standard  deviation  of  the  data 
itself,  and  multiplying  by  the  square  root  of  the  mean  thus  scaling  the  variance  to  the 
mean  of  the  digital  counts  and  enforcing  Poisson  statistics: 


dphotonsi.'E i  ?/)  ^x,y[d  counts  (x,y)\ 


d, 


counts 


[x,y) 


st  d\dcounts  (x,  y)] 


(4.1) 


Both  the  simulated  and  collected  images,  were  scaled  in  this  way  according  to  the  data 
measured  from  the  camera,  thus  converting  the  units  to  photons  for  the  reconstruction 
algorithm.  Any  read  noise  or  other  additive  white  Gaussian  noise  is  assumed  to  be 
measured  in  the  background.  These  simulated  or  experimental  images  were  then 
summed  down  each  row  to  use  the  projection  algorithm  discussed  in  section  3.1. 


f.2  Laboratory  Calibration  of  Sources 

An  experiment  was  set  up  to  mimic  the  schematic  in  figure  3.2.  Two  diodes 
were  chosen  to  provide  lOOnm  separation.  The  diodes  were  found  to  be  approximately 
560nm  and  645nm  (red  and  green)  using  a  simple  Young’s  double-slit  experiment 
setup.  In  Young’s  double-slit  experiment,  two  slits  are  used  to  form  a  diffraction 
pattern  similar  to  the  one  in  figure  4.2.  The  wavelength,  the  slit  separation  and  the 
distances  between  the  screen  and  the  slits  are  related  by: 


L  = 


A  z 
a 


(4.2) 


where  a  is  the  slit  separation,  L  is  the  distance  between  the  maxima,  z  is  the  distance 
from  the  slits  to  the  screen  and  A  is  the  center  wavelength  of  the  source  (from  [10]). 
This  equation  can  be  rearranged  to  give  the  wavelength  as  a  function  of  the  fringe 
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Figure  4.1:  Young’s  double-slit  experiment  setup 
separation  (L)  and  ratio  of  the  slit  distance  to  lens  distance  (-). 

A  =  L-  (4.3) 

z 

A  laser  with  a  wavelength  of  632. 8nm  was  used  to  determine  the  distance  ratio  - 
(hgure  4.2).  Images  were  taken  and  the  fringe  separation  L  was  calculated.  Figure 
4.3  shows  the  measurement  of  this  ratio  for  the  calibration  wavelength.  The  distance 
ratio  was  then  calculated  to  be  79.1  nm/pixel.  This  ratio  used  with  other  diodes  to 
determine  the  wavelength  center  of  each  diode.  Figure  4.4  shows  the  fringe  pattern 
from  the  higher  wavelength  diode  chosen  for  this  experiment.  Because  the  fringes  are 
saturated,  it  is  not  possible  to  find  the  center  peak.  However  the  distance  between  the 
outermost  fringe  peaks  can  clearly  be  seen.  The  distance  between  the  two  outer  peaks 
was  found,  and  the  average  period  calculated.  It  can  be  shown  from  hgure  4.5  that 
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Figure  4.2:  Fringes  for  calibration  source  at  A=632.8nm 

the  average  fringe  period  distance  L  =  8.17  pixels  and  from  equation  (4.2),  the  diode 
wavelength  is  found  to  be  A  =  645nm.  Only  the  center  wavelength  can  be  determined 
from  a  simple  double-split  experiment.  It  is  not  possible  to  determine  an  accurate 
line-width  with  this  setup.  Therefore,  the  diodes  were  chosen  with  center  wavelengths 
far  enough  apart  so  the  bandwidth  would  not  interfere  with  the  verification  of  spectral 
resolution. 

4-3  Determination  of  Lens  Dispersion 

In  this  section,  a  Focus  Aberration  Detection  (FAD)  algorithm  is  adapted  to 
determine  the  amount  of  focus  aberration  present  at  each  wavelength  using  an  algo¬ 
rithm  developed  for  a  polarimeter  using  phase  diversity  (a  defocus  aberration)  [30]. 
Since  the  diversity  element  between  this  research  and  the  polarimeter  research  is  a 
defocus  aberration,  the  FAD  algorithm  is  easily  adapted  to  the  determination  of  the 
lens-dispersion  parameter.  This  algorithm  pre-computes  the  focus  aberrations  and 
uses  a  step  similar  to  the  Richardson-Lucy  algorithm  to  deconvolve  the  focus  from 
the  object.  Finally,  the  log-likelihood  ratio  is  computed  and  a  probability  of  the  focus 
aberration  is  applied  computing  the  Maximum-a-Posteriori  probability  of  the  focus 
aberration.  In  this  implementation,  because  the  approximate  focal  length  is  known 
to  within  a  certain  range,  only  the  Maximum-Likelihood  estimation  is  needed. 
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Relative  Intensity 


Position  (pixels) 

Figure  4.3:  Calibration  fringe  distance  was  found  to  be  L  =  8  pixels 
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Relative  Intensity 


Figure  4.5:  Data  to  determine  diode  wavelength 
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The  modification  to  the  FAD  algorithm  has  4  steps.  First  start  by  generating 
a  series  of  phases  between  minimum  aberration  and  some  upper  limit  of  focus  in 
this  case  chosen  to  be  between  minimal  and  maximal  focal  lengths  observed  in  the 
lab  for  the  LEDs  above.  Secondly,  deconvolve  the  object  from  the  measured  image. 
Third,  calculate  the  log  likelihood  for  each  focus  aberration.  Finally,  find  the  maxi¬ 
mum  likelihood  which  corresponds  to  the  focus  aberration  for  the  lens  at  the  specific 
wavelength. 

A  PSF  for  defocus  is  given  by: 


Af 


hk(k,p,  A)  = 


E  A(m,  n)  exp 


7i -j(m2+n2)  (  \2Zd 


^2 


A  \N2dx 2  ( a\+/3fi)N2dx 2  J 


exp 


—‘2nj(  km+np) 
N 


Af 


E  A(m,  n)  exp 


7 rj(m2+n2)  /  A2zd _ A2zd  4 

A  \N2dx 2  (a\+/3fi)N2dx2  J 


E  exp 

x,y 


—2irj(km+np) 

N 


(4.4) 

(the  same  as  equation  (2.32)).  From  this  PSF,  the  defocus  Wd  is  varied  over  the 
distances  listed  in  table  4.1  around  the  focus  for  the  center  wavelength  Ao-  This  is 
done  for  both  wavelengths  determined  above. 
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The  log  likelihood  for  a  single  defocus  frame  is  given  by: 


LD(d )  =  ln[PD{d)\  =  EEE  dk(x,y)ln[ik(x,y)\  —  ik(x,y)  -  ln[dk(x,y)\]  (4.5) 

y  x  k 


(the  same  as  equation  (3.1)). 

The  equation  for  Wd  is  given  by: 


(4.6) 


(the  same  as  equation  (2.28)),  and  using  the  two  wavelengths  determined  above,  the 
lens  model  parameters  a  and  f3fi  can  be  computed.  The  parameters  are  listed  in  table 
4.1  given  above.  Figures  4.6  and  4.7  show  the  results  of  the  laboratory  calibration. 
Notice  that  the  PSFs  for  the  645nm  LED  do  not  visually  match  as  well,  however  the 
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Point  Spread  Function  (esimated  vs  deconvolved) 


Figure  4.6:  PSF  estimated  vs  deconvolved  for  A  =  560nm 

defocus  of  the  log-likelihood  function  was  found  to  match  closely.  This  is  due  to  the 
fact  that  although  defocus  is  the  main  aberration  present  in  the  lens,  some  of  the 
higher  order  aberrations  may  also  be  present.  The  energy  away  from  the  peak  has  a 
much  greater  overlap  and  thus  this  correlation  is  much  higher  in  the  first  (A  =  560nm) 
LED. 
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Point  Spread  Function  (estimated  vs  deconvolved) 


Figure  4.7:  PSF  estimated  vs  deconvolved  for  A  =  645nm 
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V.  Simulation  and  Experiment  Results 

The  purpose  of  this  chapter  is  to  describe  a  simulation  and  laboratory  verification  of 
the  resolution  limits  of  lens-based  chromotomagraphic  hyperspectral  sensor.  In  order 
to  accomplish  this,  a  simulation  was  performed  using  to  determine  the  effect  of  num¬ 
ber  of  defocus  planes  on  spectral  resolution.  These  simulations  were  performed  for 
multiple  wavelength  spacing  and  varying  numbers  of  defocus  planes.  In  the  previous 
chapter,  comparisons  of  various  parameters  for  the  CRLB  indicated  that  the  largest 
factor  affecting  lens-based  CTHIS  resolution  is  the  spacing  and  number  of  defocus 
planes.  Each  increase  in  the  number  of  defocus  planes  in  the  simulation  was  accom¬ 
panied  by  a  corresponding  reduction  in  integration  time,  thus  the  total  energy  in  each 
set  of  collected  images  is  the  same.  In  section  3.1,  a  projection  based  reconstructor 
was  developed  similar  to  a  previous  method  [12]  that  solves  for  the  spectrum  of  the 
scene  in  a  single  spatial  dimension  and  single  spectral  dimension.  Next,  in  section  3.2, 
the  setup  for  the  laboratory  and  the  experiment  are  described.  Finally,  a  laboratory 
experiment  performed  to  verify  a  single  instance  of  the  CRLB  and  simulation  for  the 
lens-based  CTHIS  is  described  in  section  5.4. 

5.1  Effects  of  Lens-Based  Dispersion  on  Spatial  and  Spectral  Resolution 

As  mentioned  in  section  3.1,  the  reconstruction  is  given  as  an  iterative  EM 
estimator  that  attempts  to  maximuze  the  likelihood.  The  object  estimator  began  to 
converge  around  200  iterations.  The  image  photon  bias  /3  was  chosen  to  be  10  so 
that  the  estimator  would  converge  to  a  reasonable  solution.  As  seen  in  figures  5.5 
and  5.6,  even  this  small  amount  of  bias  adds  a  great  deal  of  noise  to  the  images,  and 
results  in  a  longer  convergence  time  for  the  estimator.  Figure  5.1  gives  an  example 
of  object  reconstructions  without  and  with  bias  (using  only  1  spatial  dimension). 
The  results  of  bias  on  the  reconstruction  are  easily  seen.  The  peaks  can  be  clearly 
shown  in  the  reconstruction  with  bias,  although  there  are  many  false  peaks  and  the 
spectrum  is  blurred  across  many  more  wavelengths.  Even  though  the  of  the  bias 
reconstructed  spectrum  is  blurred,  there  are  peaks  at  the  expected  locations.  The 
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Estimated  Obje 


Figure  5.1:  Example  Reconstructions  without  (left)  and  with  (right  Bias  after  200 
iterations 

spatial  reconstruction  converges  quickly  to  the  locations  of  the  point  sources  (requiring 
only  a  few  iterations),  ffowever  the  spectrum  takes  many  more  iterations  to  show 
obvious  peaks. 

The  CRLB  is  used  as  a  metric  to  describe  the  spatial  and  spectral  resolution  of  a 
system.  Specifically,  we  use  the  standard  deviation  (the  square  root  of  the  CRLB),  and 
where  our  standard  deviation  is  lower,  the  resolution  is  considered  higher  [19] .  This  is 
equivalent  to  saying  that  when  our  estimator  standard  deviation  is  greater  than  our 
actual  point  source  separation,  it  is  impossible  to  resolve  the  object.  The  lens-based 
implementation  differs  from  other  CTHIS  implementations,  and  using  the  CRLB  is 
not  as  straightforward  as  a  spatial  resolution  criteria.  In  previous  work,  a  CRLB  was 
used  as  an  estimate  of  the  resolution  of  two  variations  of  CTHIS  sensors  [19].  The 
prism  and  grating  CTHIS  configurations  correlate  spatial  and  spectral  information 
across  the  resulting  images,  and  therefore  an  ambiguity  of  the  actual  spatial  resolution 
results.  Because  a  change  in  spatial  location  in  a  resulting  image  may  be  due  to  a 
spatial  variation  of  the  extended  object,  or  due  to  a  different  wavelength  feature,  the 
variance  of  the  spatial  separation  estimates  were  compared  with  the  actual  wavelength 
separation.  This  showed  that  the  spatial  resolution  of  the  CTHIS  sensor  for  the 
grating  and  prism  configurations  varied  with  the  spatial  separation  of  the  object. 
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Figure  5.2:  Spatial  Resolution  Estimates  vs.  Spectral  Separation 


However,  in  the  case  of  the  lens-based  CTHIS,  because  the  spectral  diversity 
varies  across  multiple  images,  we  do  not  expect  the  spatial  resolution  to  be  different 
than  a  diffraction-limited  optical  system.  Figure  5.2  shows  the  standard  deviation  and 
simulation  results  across  varying  values  of  both  spectral  separation  and  the  dispersion 
parameter  a.  The  lower  bound  obviously  bounds  the  performance  of  the  simulation, 
and  matches  the  general  trend.  The  CRLB  is  relatively  flat  with  a  changing  spectral 
separation.  Changing  the  object’s  spatial  separation  also  does  not  appear  to  affect  the 
standard  deviation.  Almost  all  of  the  simulated  points  are  below  the  Rayleigh  limit 
(the  minimum  Rayleigh  limit  for  this  system  is  2.5  pixels  at  A  =  410nm).  The  CRLB 
is  only  a  bound  on  the  best  estimator  performance,  and  is  not  necessarily  achievable. 
When  the  calculated  standard  deviation  is  below  the  Rayleigh  limit,  the  defocused 
images  must  be  providing  more  information  than  a  typical  imaging  system.  There 
may  be  special  object  cases  where,  for  this  configuration,  beyond  diffraction-limited 
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resolution  is  possible.  For  instance,  the  object  under  consideration  with  two  points 
spaced  closely,  but  distant  spectrally  obviously  benefits  from  the  increased  information 
provided  by  the  extra  frames  (which  happens  to  be  the  object  we  are  estimating).  In 
this  case,  it  would  be  obvious  there  are  two  points  present  due  to  the  fact  that  there 
will  be  almost  mutually  exclusive  zones  where  one  point  is  in  focus,  and  the  other 
point  is  out  of  focus  (see  figure  5.5  for  an  example).  Simply  observing  that  there  are 
two  zones  where  points  are  in  focus  renders  the  conclusion  that  there  are  two  points. 
This  problem  is  much  more  difficult  when  the  object  has  many  closely  spaced  points 
at  various  wavelengths  because  there  is  no  obvious  distinction  between  the  defocused 
images.  Although  the  trend  of  the  simulation  is  not  flat,  there  is  no  obvious  trend 
with  the  changing  spectral  or  spatial  separation.  As  a  result,  we  conclude  that  for 
an  average  object  that  the  best  spatial  separation  we  can  expect  from  a  lens-based 
CTHIS  would  be  around  the  same  as  the  diffraction-limited  optical  system. 

The  spectral  resolution  is  of  significant  interest.  Figure  5.3  shows  the  standard 
deviation  of  spectral  resolution  estimates  against  the  spectral  dispersion  parameter 
a.  There  is  a  trend  in  both  the  simulation  and  the  CRLB.  As  a  varies,  there  is  a 
minima  where  both  the  CRLB  and  the  simulation  give  the  lowest  <taa  before  and 
after  which,  the  standard  deviation  increases.  This  minima  occurs  for  each  spectral 
separation  A\,  as  it  increases  this  minima  requires  a  stronger  dispersion  (a  higher 
a).  When  A  a  =  lOOrtrn,  the  minima  occurs  for  a  =  80,000,  when  Aa  =  GOrtrn,  the 
minima  occurs  at  or  near  a  =  200,  000,  and  beyond  this  (for  Aa  <  40nm),  the  minima 
occurs  for  a  beyond  the  simulation  and  CRLB  parameters.  This  indicates  that  for  a 
desired  spectral  separation  (object),  there  is  a  particular  chromatic  aberration  that 
will  resolve  the  object,  but  beyond  which  the  resolution  will  get  worse.  In  other  words, 
for  a  particular  lens  with  a  known  chromatic  aberration,  there  is  a  “sweet  spot”  for 
its  spectral  resolution  showing  the  best  performance  we  can  expect  for  the  CTHIS 
system.  That  minima  is  the  limit  of  performance  for  the  given  optical  setup,  given 
the  particular  object. 
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Figure  5.3:  Spectral  Resolution  Estimates  vs.  Spectral  Dispersion 

The  number  of  defocused  planes  and  the  distances  (zj)  also  affects  the  standard 
deviation  of  the  estimator,  and  these  results  are  discussed  in  5.4,  however  for  this 
initial  study  it  was  fixed  at  20  (planes)  and  varied  ±lcm.  Increasing  the  number  of 
defocused  planes  would  require  a  precision  stage  with  which  to  move  the  lens  with 
respect  to  the  focal  plane,  but  it  was  determined  that  it  is  relatively  simple  to  move 
the  stage  in  approximately  1mm  increments.  We  would  expect  the  a  parameter  to 
dictate  an  appropriate  number  of  defocus  planes,  and  in  turn,  the  effects  on  spectral 
resolution.  These  results  mentioned  in  [18],  became  the  basis  for  the  results  reported 
in  section  5.4. 

The  Cramer-Rao  lower  bound  was  used  as  a  method  for  determining  the  ex¬ 
pected  resolution  performance  of  a  lens-based  hyperspectral  imaging  sensor  based  on 
the  lens-dispersion  parameter.  This  bound  was  verified  to  be  a  lower  bound  on  the 
expected  spatial  and  spectral  resolution  of  the  CTHIS  using  simulation.  Due  to  the 
fact  that  the  spectral  variation  is  spread  across  multiple  images,  and  not  spatially 
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distance  to  focal  plane  -  zd 


Figure  5.4:  Simplified  Optical  Setup 


as  in  the  grating  and  prism  based  CTHIS,  the  spatial  resolution  was  determined  to 
be  similar  to  that  of  a  diffraction-limited  imaging  system.  The  spectral  resolution  is 
heavily  dependent  on  the  spectral  dispersion  of  the  lens  in  the  system.  It  was  de¬ 
termined  that  each  value  of  a  yields  a  particular  “best”  point  of  spectral  resolution. 
However  it  was  also  determined  that  the  number  of  defocus  planes  was  another  pa¬ 
rameter  that  affects  spectral  resolution  and  was  not  captured  in  this  particular  use  of 
the  CRLB. 

5.2  Numerical  CRLB  Computation 

The  CRLB  given  in  section  2.4  was  computed  for  the  properties  given  in  table 
5.1.  Then,  the  resulting  images  had  a  bias  j3  added  and  noise  was  then  applied 
according  to  Poisson  statistics.  The  bias  simulates  unrelated  points  in  a  scene,  and 
also  arbitrarily  increases  the  signal  to  noise  ratio  (SNR)  of  the  images.  The  bias  was 
also  added  to  the  CRLB,  and  then  calculated  to  make  the  bound  more  realistic.  The 
final  image  then  becomes: 


ik(x,  y )  =  oihk{x  -x0,y-  y0,  A0)  +  o2hk(x  -  (. x0  +  A x),y  -  yo,  A0  +  Aa)  +  /3  (5.1) 

and  is  used  in  equation  (2.15).  Because  /5  is  a  constant,  it  only  modifies  ik(x,y )  and 
does  not  contribute  to  the  derivatives  for  the  Fisher  information  matrix  elements. 
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Figure  5.5: 
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Simulated  Defocused  Images  with  Poisson  Noise 
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Examples  of  images  with  Poisson  noise  are  given  in  figures  5.5  and  5.6  without  and 
with  a  bias  respectively. 


Parameter  Name 

Value 

A 

410nm  —  650nm 

Ao 

530 nm 

f@  A0 

200 mm 

Pfi 

158mm 

190-2  lOmrn 

#  of  defocused  planes 

3-40 

pixel  size  (dx) 

8/im 

#  of  pixels 

64 

Source  Intensity  (ol,o2) 

300  photons 

Image  Bias  /3 

10  photons 

Aperture  Radius 

2.1mm 

a 

12  x  104 

A, 

1  —  10  pixels 

Aa 

10  —  lOOnm 

Table  5.1:  CRLB  Parameters 


The  wavelengths  were  chosen  over  the  visible  spectrum  because  many  lenses 
with  severe  chromatic  aberration  are  available,  and  a  wide  variety  of  inexpensive 
detectors  are  available.  The  parameters  for  pixel  size  ( dx ),  defocused  distances  (z^), 
and  aperture  diameter  were  chosen  because  they  reflect  a  realistic  optical  system 
setup.  The  aperture  diameter  specifically  samples  at  twice  the  Nyquist  rate  for  the 
Rayleigh  criteria  of  the  pixel  size.  The  spectral  dispersion  parameter  a  was  chosen  to 
correspond  to  a  real  lens  used  in  a  laboratory  setup  (see  chapter  V).  The  maximum 
number  of  defocused  planes  k  was  fixed  at  40  because  it  corresponds  to  a  physical 
distance  of  approximately  0.5mm.  The  number  of  defocus  planes  were  varied  between 
3  and  40  and  varying  starting  positions  were  chosen  so  as  to  reduce  the  dependance 
of  the  variance  on  the  start  position,  and  to  characterize  the  resolution  dependance 
simply  on  the  number  of  frames.  The  point  source  intensities  were  chosen  for  a  very 
low  signal  setup,  and  the  image  bias  was  chosen  to  be  much  less  severe  because  it 
was  determined  that  any  higher  bias  noise  significantly  degrades  the  performance  of 
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Figure  5.6: 


Simulated  Defocused  Images  with  Poisson  Noise  and  Bias 
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an  estimator.  The  k,  Ax  and  Aa  parameters  were  then  varied  in  the  simulation  over 
the  specified  ranges  in  table  5.1. 

After  the  CRLB  was  computed  for  the  individual  images,  the  summation  over  k 
was  taken  last.  Recalling  that  equation  (2.15)  is  an  element  of  the  Fisher  information 
matrix  for  a  series  of  images,  if  cq  and  o2  are  held  constant  for  each  image,  then  the 
amount  of  light  collected  will  be  greater  the  more  frames  are  collected.  To  scale  the 
photon  count  for  the  number  images  collected,  equation  (2.15)  was  modified  to: 

p  _  1  y^  Y^  Y~^  1  dik(x,y)  dik(x,y)  ,  . 

where  K  is  the  total  number  of  frames  collected.  This  normalizes  the  number  of 
photons  across  the  all  the  images  computed  to  be  constant  regardless  of  the  number 
of  frames  collected.  This  ensures  that  collecting  more  images  does  not  result  in  a  lower 
variance  simply  because  more  photons  were  collected.  This  is  the  same  as  scaling  the 
photons  in  each  object  to  account  for  camera  integration  time. 


5.3  CRLB  Results 

As  mentioned  in  section  2.1,  the  CRLB  is  used  as  a  metric  to  describe  the  spatial 
and  spectral  resolution  of  a  system.  Specifically,  we  use  the  standard  deviation  (the 
square  root  of  the  CRLB),  and  where  our  standard  deviation  is  lower,  the  resolution  is 
considered  higher  [19].  This  is  equivalent  to  saying  that  when  our  estimator  standard 
deviation  is  greater  than  our  actual  point  source  separation,  it  is  impossible  to  resolve 
the  object.  The  standard  deviations  from  the  CRLB  a\x  and  <taa  were  then  compared 
with  Ax  and  A  a  used  to  generate  them,  and  to  the  number  of  frames  for  varying  CRLB 
calculations.  Specifically  using  the  criteria  in  equation  (2.1),  we  have 


2<?AX 

Aa 

2°aa 


1 

1 


(5.3) 
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The  number  of  frames  required  to  resolve  a  spectral  feature  were  calculated.  For 
example  in  figure  5.7,  the  line  for  A  a  =  lOrtm  never  is  never  below  20 nm  and  is 
therefore  not  resolvable  using  up  to  40  frames  to  estimate  the  scene.  It  may  be 
that  it  is  resolvable  with  greater  than  40  frames,  however  for  this  research,  the  40 
frames  was  the  maximum  number  considered.  Only  values  for  A  a  =  10  —  40nm  are 
shown  because  for  Aa  >  50 nm,  the  lines  are  too  closely  spaced  to  differentiate  them 
significantly  in  a  graph.  For  spatial  resolution,  figure  5.8  shows  an  example  of  the 


Figure  5.7:  Spectral  Standard  Deviation  versus  the  number  of  frames.  Notice  that 
the  Aa  =  lOmn  line  is  not  resolvable. 

standard  deviation  <taa  as  a  function  of  the  number  of  frames.  In  this  graph,  there 
is  also  an  unresolvable  line.  Notice  how  for  Ax  =  1,  the  standard  deviation  is  always 
greater  than  2.  This  shows  that  for  Ax  =  1  with  40  or  fewer  frames,  two  closely- 
spaced  points  are  unlikely  to  ever  be  resolvable.  These  calculations  were  made  for 
Ax  =  1  —  10  pixels  and  Aa  =  10  —  lOOnm.  The  point  where  the  ratio  in  equation 
(2.1)  is  closest  to  one,  was  calculated  and  the  number  of  frames  required  for  resolving 
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a  vs  number  of  frames,  A  =30nm 

A  A. 

X 


Figure  5.8:  Spatial  Standard  Deviation  computed  from  the  CRLB  versus  the  num¬ 
ber  of  frames.  Notice  Ax  =  1  pixel  is  not  resolvable. 

the  two  points  was  recorded.  Figure  5.9  shows  this  calculation  for  spatial  resolution. 
Only  A  a  =  10  —  40nm  points  are  shown  to  illustrate  the  trend.  This  figure  shows  that 
only  a  few  frames  more  or  less  are  needed  for  even  a  large  difference  in  Aa-  As  one 
would  expect,  this  indicates  that  the  spatial  resolution  of  a  CTHIS  system  behaves 
similarly  to  a  diffraction-limited  optical  system  and  does  not  vary  significantly  when 
points  are  spectrally  far  apart.  When  both  Ax  and  A  a  are  close  together,  many  more 
frames  are  needed  to  resolve  the  points.  When  the  points  get  spectrally  farther  apart 
(increasing  Aa),  only  a  few  frames  difference  is  required.  There  may  also  be  special 
cases  where  beyond  diffraction-limited  resolution  is  possible.  If  two  points  are  spaced 
closely  (Ax  small),  but  distant  spectrally,  then  there  is  obviously  a  benefit  from  the 
increased  information  provided  by  the  extra  frames.  It  should  be  obvious  that  two 
points  are  present  in  the  scene  if  there  are  two  mutually  exclusive  zones  where  one 


75 


#  of  frames  vs  A 


Figure  5.9:  Number  of  frames  required  to  resolve  various  spatial  separations  us¬ 

ing  the  CRLB  and  standard  deviation-based  resolution  metric.  Only  a  few  frames 
difference  are  needed  for  various  A^  values. 

point  is  in  focus,  and  the  other  point  is  out  of  focus.  The  spatial  resolution  is  improved 
by  knowing  the  spectral  characteristics  of  the  lens.  This  simple  object  however  is  not 
as  common,  and  beyond  diffraction-limited  resolution  should  not  be  normal  given  an 
average  object.  Because  of  this,  we  conclude  that  for  an  average  spectral  scene,  the 
spatial  resolution  of  the  CTHIS  should  not  be  significantly  different  from  a  diffraction- 
limited  case. 

For  spectral  resolution,  figure  5.10  shows  the  number  of  frames  required  to 
resolve  2  points  spectrally.  For  Ax  =  2  —  4  pixels,  fewer  frames  are  needed  to  resolve 
the  points  when  A^  is  small.  As  A^  increases,  fewer  frames  are  needed  to  see  the 
spectrum  of  the  object.  Notice  that  for  Ax  =  5  —  7,  only  Ax  =  7  appears  to  be  shown. 
That  is  because  for  the  CRLB,  the  profile  for  Ax  >  4  is  exactly  the  same.  In  other 
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words,  when  the  points  are  much  farther  apart  spatially,  it  is  much  easier  to  resolve 
the  points  spectrally.  For  this  case,  beyond  Ax  >  4  pixels,  there  is  no  difference  in 
the  number  of  frames  required  to  determine  the  spectrum  of  those  points.  When  Ax 


#  of  frames  vs  A. 


Figure  5.10:  Number  of  frames  required  to  resolve  various  spectral  separations. 

is  small,  many  more  frames  are  needed  to  determine  the  spectrum.  Finally,  the  closer 
the  spectrum  is  spaced,  the  more  frames  are  needed.  As  a  result,  it  seems  obvious 
that  a  coarse  resolution  spectrometer  can  be  made  using  a  very  simple  standard 
optical  setup,  for  only  a  slightly  higher  cost  of  data  collection  and  computation.  The 
chromatic  aberration  of  a  lens  is  something  that  is  present  in  many  optical  systems, 
and  the  possibility  of  collecting  more  information  for  only  a  slightly  higher  complexity 
is  very  appealing. 
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5-4  Spectral  Resolution  of  Lens- Based  CTHIS 

In  the  previous  chapter  a  resolution  criteria  based  on  the  minimum  estimator 
variance  was  presented.  The  CRLB  was  given  for  a  low-signal  scenario.  This  was  done 
because  initial  measurements  indicated  that  the  signal  to  noise  ratio  was  low.  After 
the  data  was  measured  and  converted  to  photons,  it  was  found  that  the  signal  to  noise 
ratio  was  very  high.  Also,  after  simulating  the  reconstruction  method  developed  above 
over  500  different  noise  realizations,  it  was  found  that  reconstructions  did  not  vary 
within  10  orders  of  magnitude  due  to  noise,  and  another  resolution  criteria  needed 
to  be  developed.  This  section  discusses  the  new  resolution  criteria,  and  the  results  of 
both  the  simulation  and  the  experiment. 

5-4- 1  Resolution  Criteria.  After  noisy  images  were  simulated,  the  recon¬ 
structor  developed  in  section  3.1  was  used  to  reconstruct  the  object  for  varying  lens 
positions.  Multiple  noisy  simulations  were  used  however,  the  Poisson  noise  did  not 
change  the  resulting  reconstruction  for  varying  noise  samples.  This  is  due  to  the 
fact  that  there  is  a  very  large  signal  to  noise  ratio.  Therefore,  only  1  realization  of 
noise  was  used  per  frame,  or  simply  stated,  the  simulation  was  not  repeated  multiple 
times  per  frame.  Because  of  the  lack  of  variation,  a  variance  across  multiple  random 
realizations  frames  could  not  be  determined. 

In  the  experiment,  positions  were  taken  between  the  420.4mm  and  469.2mm 
from  the  camera.  These  were  taken  to  give  equal  distances  between  the  focal  lengths 
for  the  550nm  and  the  645nm  and  equal  distances  on  either  side  of  the  focal  length.  A 
total  of  20  positions  were  taken  representing  2.5mm  which  was  the  smallest  distance 
that  could  be  taken  with  the  measurement  device.  In  simulation,  the  start  position 
was  taken  between  420.4mm  and  422mm  and  subdivided  40  times,  to  average  for  any 
resolution  criteria  no  matter  how  many  frames  were  used  from  3-20.  A  minimum  of 
3  frames  were  taken  so  there  was  always  1  frame  on  either  side  of  the  focal  point  of 
the  two  LED  wavelengths.  Although  more  wavelengths  were  used  in  simulation,  the 
experiment  only  had  2  wavelengths  represented. 
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Spectral  Energy  parameter  yx  vs.  wavelength,  3  frames,  Aj=30nm 


Wavelength  X  -  m 


Figure  5.11:  Simulated  reconstruction  for  AA=30nm  using  3  frames 

After  the  images  were  simulated  for  this  criteria,  the  objects  were  reconstructed. 
Since  we  are  primarily  concerned  with  spectral  resolution  for  this  study,  only  the  spec¬ 
tral  energy  parameter  7a  was  actually  compared  with  the  actual  spectral  separation 
Aa-  Multiple  iterations  were  attempted,  and  it  was  found  that  2000  iterations  worked 
best  for  both  the  simulated  and  experimental  data.  On  average  more  iterations  did 
not  help  the  reconstruction  further  and  in  some  cases  increased  the  amount  of  noise 
present.  When  examined,  reconstructions  using  only  a  few  frames  had  multiple  false 
peaks.  For  this  simulation  and  experiment  it  is  known  that  there  are  only  2  points. 
Figure  5.11  shows  an  example  of  multiple  peaks  for  3  frames.  The  primary  peak  should 
be  at  A  =  550nm,  and  another  peak  is  expected  at  A  =  580nm.  Notice  the  multiple 
peaks  around  the  A  =  550nm,  and  that  the  peak  at  A  =  580nm  is  almost  twice  as  large 
as  the  first  peak,  this  is  because  the  energy  for  the  first  point  is  spread  across  multiple 
wavelengths  and  more  information  is  needed  for  an  accurate  reconstruction.  Notice 
how  when  more  frames  are  added,  as  in  figure  5.12,  the  reconstruction  algorithm  is 
able  to  find  the  wavelengths  of  the  two  peaks,  even  down  to  30nm  in  simulation.  The 
peak  at  A  =  550nm  has  almost  twice  the  number  of  photons  as  the  second  peak  at 
580nm.  So  for  this  analysis  the  number  of  peaks  is  used  as  a  resolution  criteria  and 
when  two  peaks  are  detectable.  When  only  two  peaks  are  detected,  then  the  image 
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Spectral  Energy  parameter  y^  vs.  wavelength,  7  frames,  A^=30nm 


Figure  5.12:  Simulated  reconstruction  for  A_\=30nm  using  7  frames 

is  considered  spectrally  resolved.  This  resolution  criteria  is  in  keeping  with  the  FAD 
algorithm  discussed  earlier,  as  well  as  other  iterative  maximum-likelihood  estimators 
in  the  literature  where  the  values  for  which  the  algorithm  converges  are  considered 
the  correct  parameter  for  estimation.  In  the  case  of  the  FAD  algorithm,  the  amount 
of  focus  error  is  considered  a  correct  estimate  when  the  estimation  converges.  Fig¬ 
ure  5.11  represents  an  unresolved  case  and  figure  5.12  represents  a  resolved  case.  A 
modified  first  derivative  test  was  used  to  detect  peaks.  The  first  derivatives  were  taken 
and  the  zero  crossings  were  determined  to  find  the  critical  points.  The  zero  crossing 
location  were  compared  to  determine  if  they  were  an  inflection  point,  a  minima  or  a 
maxima.  If  a  maxima  occurs  at  the  critical  point,  the  peaks  were  kept  if  they  were 
within  2  orders  of  magnitude  of  the  known  intensity  for  the  second  peak  (in  this  case 
at  580rtrn) .  Specifically  the  threshold  is  at  5  x  1011  for  this  case. 

5-4-2  Simulation  and  Laboratory  Results.  Images  were  simulated  and  col¬ 
lected  as  described  above  for  each  set  of  different  frames.  The  total  number  of  frames 
were  taken  for  each  start  position,  then  frames  were  successively  added  in  and  each  set 
then  was  reconstructed  using  the  algorithm  developed  in  section  3.1,  and  the  modified 
second  derivative  test  applied  to  determine  the  peaks.  The  number  of  peaks  over  each 
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Simulated/Experimental  Mean  number  of  peaks  found  A^=80-100nm 


Figure  5.13:  Mean  Number  of  peaks,  simulated  and  experimental 

start  position  was  averaged  for  the  simulation,  and  averaged  using  a  weighted  average 
for  the  experimental  data.  Figure  5.13  shows  how  the  number  of  frames  affects  the 
average  number  of  peaks  for  80 nm  to  lOOnm.  Although  the  data  was  simulated  for 
10 nm  to  lOOnm  only  80 nm  to  lOOnm  is  displayed  because  this  region  displays  the 
clearest  trend.  As  the  wavelength  increases,  the  number  of  peaks  slightly  decreases. 
In  other  words,  the  farther  away  two  points  are  in  spectrum,  the  better  the  recon¬ 
struction  is,  with  fewer  frames.  Because  the  trend  is  so  strong  with  only  the  first 
few  frames  between  3  frames  and  7  frames,  improvement  due  to  increased  A  a  is  only 
slight,  however  it  does  demonstrate  that  the  farther  apart  two  wavelengths  are  being 
reconstructed,  the  easier  it  is  to  resolve  them.  Notice  also  that  after  7  or  8  frames, 
almost  all  of  the  wavelengths  are  resolvable.  At  this  point  we  have  demonstrated  the 
trend  with  wavelength  on  resolvability,  we  now  seek  to  understand  how  closely  the 
measured  A  a  relates  to  the  actual  Aa.  For  the  experimental  data,  we  see  the  same 
trend:  the  farther  apart  the  points,  the  fewer  frames  that  are  needed  to  resolve  the 
data.  Notice  in  figure  5.14  of  the  experimental  data  that  the  two  points  are  clearly 
separated  unlike  figure  5.11  and  are  therefore  resolvable.  However,  the  intensities  are 
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Spectral  Energy  parameter  yk  vs.  wavelength,  3  frames,  A^=95nm 


Figure  5.14:  Experimental  reconstruction  for  A^=95nm  using  3  frames 

not  the  expected  values.  The  first  peak  is  at  550nm  as  expected,  but  the  second 
peak  is  at  620nm  which  is  only  Aa  =  70nm  when  the  actual  value  is  Aa  =  95nm 
as  determined  by  the  Young’s  experiment  described  in  section  4.2.  The  data  is  re¬ 
solvable  with  only  3  frames,  however  the  intensity  and  the  spectral  accuracy  are  still 
not  what  we  expect.  Figure  5.15  shows  for  an  example  of  the  experimental  data 
that  for  7  frames,  the  intensities  are  closer  to  their  expected  values,  especially  their 
relative  values.  However  the  spectral  accuracy  is  still  only  Aa  =  80nm  rather  than 
Aa  =  95nm  as  expected.  As  another  example  using  the  full  20  frames  reconstruction, 
the  intensities  are  exactly  where  they  should  be,  and  the  points  are  clearly  resolvable, 
however,  Aa  =  70 nm  rather  than  Aa  =  95 nm.  The  discrepancy  in  the  measured 
spectral  accuracy  with  the  calibration  vs.  the  reconstructed  data  may  be  due  to  the 
measurement  error  in  the  Young’s  calibration,  or  most  likely  due  to  the  sensitivity 
of  the  reconstructed  wavelength  bins.  Point  spread  functions  were  generated  every 
lOnm  for  the  reconstruction,  and  there  might  be  some  sensitivity  in  the  experimental 
data  to  an  optimal  wavelength  spacing  for  PSFs  compared  to  the  distance  the  lens 
is  moved  in  between  frames.  Since  the  data  was  simulated  for  multiple  wavelengths, 
the  spectral  accuracy  can  be  computed  from  the  simulated  data.  Figure  5.17  shows  a 
graph  of  the  average  spectral  accuracy  (computed  Aa)  as  a  function  of  the  number  of 
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Spectral  Energy  parameter  vs.  wavelength,  7  frames,  A^=95nm 


Figure  5.15: 


Figure  5.16: 


Wavelength  X  -  m 


Experimental  reconstruction  for  A^=95nm  using  7  frames 


Wavelength  X  -  m 


Experimental  reconstruction  for  A^=95nm  using  20  frames 
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^  o_8  Simulated  Spectral  Accuracy  Mean  A^-nm,  Ax=80-1  OOnm 


Figure  5.17:  Mean  simulated  spectral  accuracy  for  80-100nm  vs.  number  of  frames 

frames.  Notice  that  although  the  data  is  close  to  the  simulated  A^  value,  the  values 
are  around  5 nm  away  from  the  actual  value  even  after  the  data  is  resolvable.  Another 
metric  that  can  be  calculated  when  examining  spectral  accuracy  is  the  sum  squared 
error.  This  is  calculated  by  squaring  the  difference  between  Aa  and  Aa  the  estimated 
value.  Figure  5.18  shows  that  the  error  drops  off  significantly  after  6  —  7  frames.  At 
less  than  6  frames  the  greater  the  separation  Aa,  the  lower  the  error.  Figure  5.19 
that  the  CRLB  bounds  the  performance  of  the  simulated  and  experimental  results 
and  provides  a  metric  for  determining  spectral  resolution  performance,  ft  should  be 
noted  that  the  CRLB  has  more  information  than  the  reconstructor  due  to  the  fact 
that  the  object  structure  (2  points)  is  known  for  the  CRLB.  This  is  required  in  order 
to  compute  the  CRLB,  but  although  a  two-point  reconstructor  could  be  developed 
it  would  not  be  a  realistic  algorithm  for  experimental  use  as  it  could  only  detect 
two-point  objects. 
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Simulated  Spectral  Accuracy  Sum  Squared  Error  ZA  ,  A^=80-100nm 


x  10' 


Figure  5.18:  Sum  Squared  Error  simulated  spectral  accuracy  for  80-100nm  vs. 

number  of  frames 


#  of  frames  vs 


Figure  5.19:  CRLB  Plotted  with  Simulation  and  Experimental  results 
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5.5  Lens  Point  Spread  Function  Modelling  Signal-to-Noise  Ratio  as 
Compared  to  the  CRLB 

In  calculating  real  results,  it  is  necessary  to  compare  the  models  to  laboratory 
data.  Since  the  model  does  not  match  the  real  noise  also  received  in  the  image,  the 
simplified  point  spread  function  does  not  account  for  everything  from  the  real-world 
scenario.  This  section  will  discuss  the  quantification  of  this  modelling  error  and  lay 
the  foundation  for  comparing  the  Cramer- Rao  lower  bound  to  the  real  experiment. 
The  real  PSF  in  the  presence  of  modelling  error  is  given  by: 


h(x,  y)  =  (g(x,y)+e(x,  y)) 


(5.4) 


where  g(x,y)  is  the  modelled  PSF  and  the  image  resulting  from  this  noisy  PSF  (with 
an  ideal  object)  is: 

i(x,y)  =  h(x,y)  *  o(x,y)  (5.5) 

As  we  can  see  from  equation  (5.4),  the  PSF  has  modelling  error  in  it  that  is  not 
accounted  for  by  the  model.  This  error  is  then  present  in  the  received  images  used 
for  calibrating  the  PSF  of  the  lens.  Using  the  deconvolution  mentioned  above,  the 
estimated  point  spread  function  h'(x,y )  is  calculated  which  was  then  used  to  deter¬ 
mine  the  parameters  of  the  lens.  To  compare  the  model  to  the  real  PSF  we  start  by 
normalizing  the  power  in  the  PSF  to  unity: 


hd{x,  y) 


h'(x,  y)  -  E[h'(x,y)\ 
Y,xh'{x,y) 


(5.6) 


Then,  the  normalized  PSF  hd(x,y )  is  convolved  with  an  ideal  object  to  come  up  with 
an  estimate  of  the  image  with  error: 


ii(x,y)  =  hd(x,  y)  *  o(x,y) 


(5.7) 


where  o(x,  y )  is  the  same  object  used  in  the  simulation  of  the  CRLB  and  the  simulation 
and  we  call  the  image  with  error  id(x,y).  Next,  the  modelled  PSF  g(x,y )  convolved 
with  the  ideal  object  to  calculate  an  estimate  of  the  error-free  image. 


i(x,y)  =  g(x,y)to(x,y) 


(5.8) 


Finally  the  estimate  of  the  modelling  signal  to  noise  ratio  is  computed  using  the 
mean  of  the  ratio  of  the  image  divided  by  the  error.  This  is  done  for  some  small 
distance  around  the  peak  (a  square  of  10  pixels)  to  capture  signal,  and  leave  out  the 
background. 


SNR 


model 


=  Er 


i(x,y) 


(5.9) 


x'y  L \Kx>y) 

The  SNR  of  the  CRLB  is  calculated  using  the  intensities  of  the  sources  divided  by 


the  diffuse  background  intensity. 


SNRcrlb  —  Ex,y 


i(x,y)~ 

VB  _ 


(5.10) 


In  this  case,  the  SNRmodei  =  8.24  and  with  o\  =  o2  =  300  and  B  =  10,  the 
SNRcrlb  =  3.59.  This  compares  well  and  shows  that  even  though  the  intensity 
of  the  experiment  was  much  higher  than  the  signal  of  the  CRLB,  the  SNR  of  the 
modelling  error  for  the  experiment  compares  well  with  the  SNR  of  the  CRLB. 


5.6  CRLB  as  a  Metric  for  Lower  Bound  on  Spectral  Resolution 

This  chapter  has  examined  a  simulation  and  laboratory  verification  of  spectral 
resolution  and  also  examined  factors  affecting  spectral  accuracy.  The  primary  factor 
in  this  study  was  determined  to  be  the  number  of  frames.  There  appears  to  be  a 
minimum  number  of  frames  required  in  order  to  resolve  two  closely-spaced  spectral 
points,  beyond  which  more  frames  gives  little  to  no  additional  information.  The 
farther  apart  these  points,  the  fewer  the  frames  necessary  to  resolve  two  points,  and 
the  greater  the  spectral  accuracy.  A  simulation  was  set  up  to  determine  the  extent 
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that  these  factors  affect  spectral  resolution  and  spectral  accuracy.  An  experiment  was 
set  up  and  demonstrated  the  ability  of  the  presented  algorithm  to  extract  meaningful 
color  information  from  a  series  of  defocused  frames  using  chromatic  aberration.  Also, 
the  trend  for  more  frames  to  not  affect  spectral  resolution  and  accuracy  beyond  6 
frames  for  this  setup  was  determined  by  the  simulation  and  verified  by  the  experiment. 
The  CRLB  compares  well  with  the  simulated  and  experimental  results  as  shown  in 
section  5.1  where  the  CRLB  was  used  to  determine  what  effect  lens  dispersion  would 
have  on  spectral  performance,  as  well  as  in  section  5.4.2  where  the  parameter  study  of 
the  number  of  frames  was  compared  with  both  a  simulation  and  an  experiment.  While 
the  CRLB  does  not  completely  characterize  the  design  of  a  lens-based  CTHIS,  it  can 
be  useful  as  a  method  of  bounding  the  effects  of  various  parameters  on  reconstruction 
and  useful  for  finding  proper  parameters  to  design  a  sensor.  While  the  analysis  of  this 
chapter  has  focused  on  a  lab-based  setup  for  sensor  testing,  the  next  chapter  gives 
a  simple  study  for  the  lens-based  CTHIS  in  the  presence  of  the  atmosphere,  which 
has  the  potential  to  significantly  degrade  the  sensor  performance  and  makes  it  more 
difficult  to  use  these  CTHIS  sensors  in  real-world  applications. 


VI.  Blind  Deconvolution  and  Hyperspectral  reconstruction 
in  the  presence  of  Atmosphere 
6.1  Atmospheric  Results 

In  this  section,  the  algorithm  given  in  section  3.1  is  re-derived  for  a  CTHIS  in 
the  presence  of  the  effects  of  atmospheric  turbulence.  This  is  a  blind  deconvolution 
approach  which  attempts  to  determine  the  seeing  parameter  of  the  atmosphere  from 
a  series  images  collected  with  a  known  amount  of  chromatic  aberration  produced  by 
a  lens.  The  differences  for  the  assumptions  from  a  laboratory  setup  are  discussed 
in  section  6.2.  Section  6.3  discusses  how  the  Expectation-Maximization  approach 
for  estimating  a  spectral  and  spatial  scene  derived  in  section  3.1  can  be  extended  to 
jointly  estimate  the  background  and  the  spatial-spectral  object  in  the  presence  of  an 
atmospheric  Optical  Transfer  Function  (OTF)  and  also  used  to  estimate  the  seeing 
parameter.  Finally,  section  6.4  discusses  the  simulation  of  two  scenarios  showing  the 
application  of  this  algorithm. 


6.2  CTHIS  System  Design  and  Modelling 

The  model  used  in  (2.26)  was  strictly  based  on  an  assumption  of  the  focal 
length.  A  more  accurate  lens  model  for  a  lens-based  CTHIS  incorporates  the  index 
of  refraction  change  with  wavelength  n{ A).  The  wavelength  dependence  of  the  focal 
length  of  a  typical  uncorrected  thin  lens  with  equal  radii  of  curvature,  R,  for  the  front 
and  back  surfaces  is  given  by  equation  (6.1). 


/(A/) 


R 

2n{Xi)  -  2 


(6.1) 


In  this  equation  A /  is  the  wavelength  of  the  light  passing  through  the  lens  associated 
with  a  discrete  spectral  component  indexed  by  l  and  n  is  the  index  of  refraction  [11], 
As  previously,  the  focal  length  is  assumed  known  as  a  function  of  wavelength,  and  it 
is  used  to  define  the  known  point  spread  function  (PSF)  of  the  optical  system,  hopt, 


which  is  approximated  by  equation  (6.2). 


hopt(x,y,l,k ) 


N  N 

A(w,  u)e  A 

W=1  V=1 


il  (it2-  V2)All2  _3^(ux+vy)AXAu 

f(xi)Jy  '  e  x‘Zk 


(6.2) 


This  is  the  same  as  previously  assumed  in  (2.24)  where  Zy.  is  the  distance  between  the 
lens  and  the  detector  plane.  The  PSF  model  presented  in  equation  (6.2)  is  a  fairly 
good  approximation  of  the  true  function  when  the  object  being  imaged  is  far  enough 
away  from  the  system  to  be  considered  at  infinity  [11].  The  distance  Zk  from  the 
lens  to  the  plane  where  the  image  is  formed  for  the  kth  image  taken  by  the  sensor, 
A  is  the  pupil  function,  (x,y)  are  sample  locations  in  the  detector  plane,  (u,v)  are 
sample  locations  in  the  pupil  plane,  Au  is  the  sample  spacing  in  the  pnpil  plane  and 
A*  is  the  sample  spacing  in  the  detector  plane.  A  discrete  approximation  of  the  PSF 
is  needed  to  model  the  PSF  in  a  computer.  The  OTF  ( Hopt )  is  the  2-dimensional 
Fourier  transform  of  the  PSF  generated  using  equation  (6.2). 

Another  diffraction  related  effect  is  contributed  by  atmospheric  turbulence  and 
can  be  modelled  using  an  average  optical  transfer  function  model.  The  effective 
average  transfer  (Hatm)  in  equation  (6.3)  is  for  an  optical  system  viewing  an  object 
through  a  long-exposure  image  viewed  through  turbulent  atmosphere  with  a  seeing 
parameter  of  r0(Xi)  [13]. 


H, 


atm 


(6.3) 


In  this  equation  ( fx,fy )  are  coordinates  in  frequency  space  [10].  The  seeing  parameter 
is  dependent  on  which  spectral  band  is  used  in  the  model  and  is  assumed  to  vary  with 
wavelength  via  the  following  function: 

roW  =  r?a(YJ-)  (6-4) 

where  A min  is  the  minimum  wavelength  in  the  pass-band  of  the  optical  filter  and  r”wn  is 
the  seeing  parameter  corresponding  to  this  wavelength.  The  atmospheric  and  optical 
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components  of  the  transfer  function  combine  to  form  the  total  transfer  function  (Htot). 
The  total  PSF  (htot(x,  y,  l,  k,  r™m)  is  the  2-D  inverse  FFT  of  Htot. 


Htotifx,  fy ,  l,  k,  Cin)  =  Hopt(fx,  fy ,  /,  k)Hatm(fx ,  fy,  l )  (6.5) 

The  images  formed  by  the  CTHIS  sensor  are  panchromatic  images  featuring  one 
wavelength  of  light  in  focus  and  others  across  the  remaining  spectral  components  out 
of  focus  to  varying  degrees.  The  optical  system  is  modeled  as  being  linear  in  intensity 
and  wavelength  allowing  the  superposition  principle  to  be  used  in  modeling  the  signal 
vectorized  by  the  charge-coupled  device  (CCD)  array.  The  equation  for  the  sampled 
projection  vector  in  units  of  photons,  /,  is  shown  in  equation  (6.6). 

I{rn,  k)  fa  °(x’  V’  l)htot(i  -x,m-  y,  l,  k,  rfn)  (6.6) 

i  l  x  y 

=  °(^>  l)htot(m  -  y,  /,  k,  r“in)  (6.7) 

i  y 

In  this  equation,  o  is  the  intensity  of  light  in  units  of  photons  falling  on  the  detector  at 
any  point  (x,y)  predicted  by  geometric  optics.  The  variables  (i,m)  are  also  detector 
coordinates  in  units  of  samples.  An  index  l  indicates  which  spectral  component  of 
the  scene  is  under  consideration  and  k  is  an  index  referencing  the  discrete  distance 
Zk  between  the  lens  and  the  image  plane.  Note  that  equations  (6.2)  and  (6.6)  are 
approximations  due  to  modelling  both  the  electro-magnetic  spectrum  and  the  images 
at  each  spectral  component  as  discrete  portions  of  the  true  intensity  seen  on  the  CCD 
array.  Because  of  the  projection  operation  and  associated  dimensionality  reduction 
the  model  can  be  simplified  in  terms  of  the  spectral  scene  vector,  o  and  the  vector 
impulse  response  of  the  system,  htot  shown  in  equation  (6.6). 

The  data  gathered  by  this  sensor  contains  random  components  as  well  as  those 
predicted  by  radiometry  and  diffraction  theory.  The  noise  present  on  the  sensor  is 
assumed  to  be  generated  by  two  components,  Poisson  noise  and  read  noise.  The  first 
is  due  to  the  random  arrival  times  of  photons  in  the  light  itself  which  is  known  to  have 
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a  Poisson  probability  mass  function.  The  second  is  the  noise  generated  from  the  dark 
current  flowing  through  the  CCD  detectors.  The  expected  number  of  dark  counts  at 
each  detector  pixel  ( B )  is  assumed  to  be  unknown  at  the  time  of  observation  and  is 
the  mean  of  a  Poisson  random  variable  representing  the  number  of  actual  dark  counts 
measured  at  each  pixel  location.  The  total  signal  D  measured  by  the  system  is  shown 
in  equation  (6.8).  C(m )  is  a  function  that  is  binary  in  nature  and  represents  an  area 
of  the  CCD  where  light  passes  to  the  detector  array.  The  purpose  of  this  detector 
aperture  function  is  to  define  an  area  where  the  signal  exists  in  order  to  facilitate  the 
estimation  of  the  photo-detector  bias  without  photons  from  the  scene.  This  can  be 
done  in  the  laboratory  with  a  stop  or  mask  placed  before  the  detector. 

D(m,  zk )  =  C(m)I(m,  zk)  +  B  +  qQ(m,  z)  +  qb(m,  z)  (6.8) 


6.3  Scene  Reconstruction  and  Seeing  Parameter  Estimation 

In  this  section,  the  reconstructor  previously  developed  in  section  3.1  is  extended. 
The  lens-based  chromotomagraphic  hyperspectral  sensor  takes  a  series  of  images  each 
having  a  known  defocus  from  the  other  images.  Given  the  collected  image  data, 
the  originating  scene  is  unknown  so  an  estimator  the  must  be  used.  In  this  case, 
the  unknown  parameters  are  the  scene  intensity  at  each  pixel  location  and  spectral 
component,  o(y,l),  the  seeing  parameter  of  the  atmosphere  ra  and  the  dark  current 
generated  bias  signal  level  at  each  detector  pixel  B. 

The  EM  algorithm  discussed  in  section  3.1  is  applicable  to  estimating  the  pa¬ 
rameters  of  interest  with  a  bit  of  modification.  Substituting  equation  (6.8)  into  the 
joint  probability  of  all  the  data  previously  derived  in  section  3.1  (equation  (3.16))  the 
joint  probability  mass  function  is  shown  in  equation  (6.9). 


N  K 


P[D  =  dV(m,  k)\  =  ]  [  f  j 


m=  1  k= 1 


(■ C(m)I(m ,  zk)  +  B)d^m’Zk^ 
d(m,zk)\ 


( C(m)I(m,Zk)-\-B ) 


(6.9) 
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The  products  in  equation  (6.9)  are  over  the  number  of  pixels  in  the  vertical  directions 
as  well  as  the  number  of  images  gathered  from  the  system,  K .  Maximizing  this  joint 
probability  as  a  function  of  the  scene,  o,  the  pixel  bias,  B  and  the  seeing  parameter  r0, 
is  challenging.  However,  using  the  EM  algorithm  (figure  3.1),  the  previous  derivations 
can  be  modified  to  be  used  in  this  situation.  Previously  the  background  was  only 
computed  or  directly  measured  as  a  constant  rather  than  estimated  each  iteration.  In 
this  setup,  the  background  was  jointly  estimated  with  the  object.  To  incorporate  a 
background  estimate  B  into  each  iteration,  we  assume  the  complete  data  are  random 
variables,  D\  and  D2  whose  means  are  given  by: 

E[Di(m,  y,  l,  k)]  =  C(m)o(y,l)htot(m  -  y,l,k,  r0mm)  (6.10) 

and 

E[D2(m,k)]  =  B  (6.11) 

where  D\  is  using  the  total  transfer  function  and  includes  the  atmospheric  seeing 
parameter  r0.  The  choice  to  define  two  distinct  sets  of  complete  data  is  not  by 
any  means  the  only  one,  but  has  been  shown  to  be  advantageous  in  other  imaging 
problems.  As  previously,  the  complete  data  are  related  to  the  incomplete  data  via 
the  following  transformation. 

N  L 

D(m,  k)  =  EE  D1(m,y,l,k)  +  D2(m,k)  (6.12) 

y= i  i=i 

This  relationship  is  statistically  consistent  if  the  complete  data  sets  are  chosen  to  be 
independent  Poisson  random  variables  since  the  sum  of  Poisson  random  variables  is 
Poisson  as  well.  Also  adding  together  Poisson  random  variables  produces  a  Poisson 
whose  mean  is  equal  to  the  means  of  the  random  variables  being  added  together.  With 
the  statistical  model  for  the  complete  data  in  hand,  and  using  the  techniques  previ¬ 
ously  discussed  in  section  3.1,  the  iterative  update  solution  for  the  spatial-spectral 
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object  o(y0,l0 )  is  given  by: 


new,  ,  )  =  old,  ,  \srsr _ d(m,  k)htot(m  —  y0,  k,  l0,  r™in) _ 

(/oWK  fc)  +  E  E  C(m)htot(m  -  Vo,  k,  la,  r“in) ' 

km 

(6.13) 

Similar  in  form,  the  joint  estimation  for  the  background  is  given  by  the  following 
update  equation: 


j^new 


K  M 

k=  1  m=l 


d(m,  k ) 

K(C(m)Iold(m,  k )  + 


(6.14) 


This  solution  has  the  advantageous  property  that  the  estimated  bias  level  is  always 
positive.  The  update  equations  in  (6.13)  and  (6.14)  are  repeated  every  iteration  until 
a  termination  criteria  is  satisfied.  In  the  case  where  photon  counting  noise  is  the 
dominant  source  of  error,  the  following  criteria  can  be  used  to  stop  the  iteration: 


Y  k )  “  C(m)Inew(m,  k)  -  Bnewf  <  Y  Y  k )  (6‘15) 

in  k  in  k 

This  criterion  is  derived  from  the  relationship  between  the  mean  and  variance  of  a 
Poisson  random  variable.  The  left  side  of  the  equation  is  related  to  the  variance  of 
the  estimated  noise  and  the  left  side  is  related  to  the  sample  mean  in  that  both  sides 
are  the  non-normalized  versions  of  the  variance  and  mean  respectively.  The  iterations 
continue  until  either  this  criterion  is  met  or  a  maximum  number  of  iterations  are 
reached. 

Thus  far  no  mention  has  been  made  of  how  to  estimate  the  seeing  parameter,  rG, 
from  the  imagery.  The  strategy  for  estimating  the  seeing  parameter  is  to  execute  the 
EM  algorithm  for  a  number  of  different  values  of  rQ  and  choose  the  smallest  value  for 
which  the  convergence  criteria  is  satisfied.  The  hypothesis  for  using  multiple  values 
of  ra  is  that  if  the  seeing  parameter  is  too  low,  the  estimated  noise  will  be  much  larger 
than  the  shot  noise  component.  The  estimator  will  then  not  be  able  to  match  the 
data  with  a  choice  of  spectral  scene  in  cases  where  the  seeing  parameter  is  too  low, 
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Table  6.1:  System  parameters  chosen  for  the  simulation 


Band-Pass  Filter  Transmission 

none 

1  or  0 

Band-Pass  Filter  Bandwidth 

nm 

500-690  (10  nm  steps) 

CD  array  Array  Size 

pixels 

128  X  128 

CCD  array  Pixel  Pitch 

mm 

5.0  (center  to  center) 

Focusing  Optics  Diameter 

Meters 

0.5  (across  the  optic) 

Focusing  Optics  Focal  Length 

Meters 

10  (for  500  nm  light) 

Focusing  Optics  Focal  Length 

Meters 

10.01  (for  700  nm  light) 

due  to  the  limited  bandwidth  of  the  impulse  response.  The  mismatch  between  system 
bandwidth  and  data  bandwidth  produces  excess  noise  which  prevents  the  algorithm 
from  converging.  Convergence  is  possible  for  larger  seeing  parameter  values.  However 
the  larger  the  seeing  parameter  is,  the  lower  the  bandwidth  of  the  reconstructed 
spectral  scene  o.  Since  we  desire  the  sharpest  spectral  scene  reconstruction  while  still 
obtaining  algorithm  convergence,  the  minimum  seeing  parameter  value  that  achieves 
the  criterion  in  equation  (6.15)  is  selected. 


6.4  Algorithm  Performance 

The  performance  of  the  proposed  algorithm  is  tested  using  simulated  binary 
source  data  featuring  different  separations  of  the  sources  in  both  wavelength  and 
space.  Different  signal  levels  are  explored  as  well  as  different  levels  of  background 
radiation.  In  this  way  the  signal  to  noise  ratio  of  the  data  can  be  controlled  as  different 
levels  of  spectral  and  spatial  resolution  are  investigated.  The  simulated  system  used 
to  generate  the  data  is  described  in  subsection  6.4.1  of  this  chapter.  The  results 
obtained  from  testing  the  algorithm  on  the  simulated  data  are  shown  in  subsection 
6.4.2. 


6-4-1  Closely  Spaced  Sources  Separated  in  Wavelength  by  100  nm.  A  specific 
set  of  system  parameters  are  chosen  for  simulating  realistic  CTHIS  data.  The  partic¬ 
ular  system  parameters  are  consistent  with  those  of  a  small  telescope  imaging  system. 
Table  6.1  contains  the  parameters  of  the  electro-optical  system.  The  simulated  system 
is  designed  to  take  20  vectors  sampled  at  regular  intervals  as  the  distance  between  the 
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X  -  p.METERS 

Figure  6.1:  2  optical  sources  in  position  (Vertical)  and  wavelength  (Horizontal) 

lens  and  the  CCD  array  changes  from  10  meters  to  10.01  meters.  The  targets  viewed 
by  the  system  are  mono-chromatic  sources  separated  by  varying  degrees  in  space  and 
wavelength.  The  first  experiment  carried  out  in  this  study  involves  the  two  sources 
placed  on  top  of  one  another,  (no  physical  separation)  but  separated  in  wavelength 
by  100  nanometers.  The  sources  each  provide  10000  photo-electrons  to  the  imaging 
system  for  each  frame  taken  with  100  photo-electrons  of  dark  current  being  read  out 
at  each  detector  pixel.  The  seeing  parameter  is  chosen  to  be  15  centimeters.  Figure 

6.1  shows  an  image  of  the  sources  as  a  function  of  wavelength  and  position.  Figure 

6.2  shows  simulated  frames  of  spectral  vectors  as  a  function  of  lens  position  away 
from  the  focal  length  of  10  m.  Figure  6.3  shows  the  reconstructed  spectral  image. 
The  algorithm  identified  the  seeing  parameter  as  being  equal  to  15  centimeters  when 
searching  on  a  range  from  10  to  20  cm. 

The  experiment  was  repeated  for  a  lower  value  of  signal  photons.  In  this  second 
case  the  photon  level  of  the  sources  was  dropped  to  1000  photons.  The  raw  projection 
data  is  shown  in  Figure  6.5. 
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LENS  POSITION  -  METERS  x  10'3 

Figure  6.2:  Simulation  of  2  optical  sources  as  viewed  through  a  lens-based  chromo- 
tomographic  imager.  The  vector  projection  readout  is  shown  in  the  columns  of  the 
image.  Each  column  corresponds  to  a  frame  of  data  taken  at  a  position  marked  on 
the  horizontal  axis  that  represents  a  deviation  from  the  10  meter  position  from  the 
primary  focusing  optic. 

The  reconstructed  spectral  projection  estimates  are  shown  in  figure  6.5.  Al¬ 
though  the  estimated  spectral  projections  are  fairly  accurate  (but  demonstrate  some 
additional  spectral  width),  the  estimated  seeing  parameter  was  15  cm.  The  error  in 
the  spectral  projection  estimates  due  to  the  lower  signal  to  noise  ratio  did  not  affect 
the  algorithm’s  ability  to  calculate  the  correct  value  for  the  seeing  parameter. 

6-4-2  Sources  Separated  Only  in  Wavelength  by  40  nm.  The  sources  each 
provide  10000  photo-electrons  to  the  imaging  system  for  each  frame  taken  with  100 
photo-electrons  of  dark  current  being  read  out  at  each  detector  pixel.  The  seeing 
parameter  is  chosen  to  be  15  centimeters.  Figure  6.1  shows  an  image  of  the  sources  as 
a  function  of  wavelength  and  position.  Figure  6.6  shows  simulated  frames  of  spectral 
data  at  the  distances  where  the  two  sources  would  be  in  focus.  Figure  6.7  shows  the 
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Figure  6.3:  Reconstructed  spectral  image  of  the  two  optical  sources  from  the  data 
presented  in  Figure  6.2.  The  seeing  parameter  was  estimated  to  be  15  cm  and  the 
true  value  was  15  cm. 

reconstructed  spectral  image.  The  algorithm  identified  the  seeing  parameter  as  being 
equal  to  15  centimeters  when  searching  on  a  range  from  10  to  20  cm. 

6.5  Discussion  of  Blind  Deconvolution  of  Hyperspectral  Data  in  the 
presence  of  Atmosphere 

The  proposed  algorithm  for  reconstructing  spectral  projections  from  chromoto- 
mographic  vector  projection  data  while  simultaneously  estimating  the  seeing  param¬ 
eter  through  which  the  image  data  is  gathered  is  demonstrated  to  work  at  signal  to 
noise  ratios  between  30  and  100.  The  algorithm  is  presumed  to  work  properly  for 
higher  SNR  conditions  but  would  require  more  time  to  converge  thus  making  that 
study  more  time  consuming.  Further  trials  need  to  be  conducted  to  determine  the 
range  of  signal  to  noise  ratios  and  achievable  spectral  resolutions  over  which  the  algo¬ 
rithm  will  perform  well.  Also,  experiments  with  measured  data  should  be  conducted 
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Figure  6.4:  Raw  projection  data  for  the  case  where  both  sources  provide  1000 

photons  during  the  measurement  time  and  the  background  level  is  set  at  100  photons. 
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Figure  6.5:  Reconstructed  spectral  projections  for  the  low  SNR  case.  The  spectral 
features  are  broadened  but  still  distinguishable. 
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Figure  6.6:  Raw  projection  data  of  two  objects  4  pixels  apart  and  separated  by  40 
nanometers  in  wavelength. 


to  demonstrate  the  utility  of  the  algorithm  in  the  presence  of  modeling  error  and 
other  unaccounted  for  but  unavoidable  effects. 
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Figure  6.7:  Reconstructed  spectral  projections  of  features  separated  by  40  nanome¬ 
ters. 
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VII.  Conclusions 


This  research  has  examined  several  factors  affecting  resolution  of  a  lens-based  chro- 
motomographic  hyperspectral  imaging  sensor  in  the  presence  of  noise.  A  lower  bound 
on  resolution  was  developed  to  use  as  an  engineering  design  tool  to  determine  how 
changes  in  sensor  design  affected  resolution.  A  reconstructor  was  derived  to  account 
for  background  noise  which  accounted  for  a  realistic  laboratory  setup.  Finally,  a  simple 
experiment  was  performed  and  used  with  the  reconstructor  to  verify  the  lower-bound 
and  its  usefulness  as  a  metric.  The  effects  of  the  amount  of  lens  chromatic  aberration 
was  studied  and  determined  to  be  a  factor,  but  the  most  significant  factor  affecting 
spectral  resolution  of  a  lens-based  CTHIS  was  found  to  be  the  number  of  defocus 
planes.  This  number  of  defocus  planes  compares  well  with  simulation  and  the  exper¬ 
iment.  Also,  a  modified  blind  reconstructor  was  derived  in  order  to  account  for  the 
presence  of  atmospheric  turbulence.  This  blind  reconstructor  estimates  jointly  the 
background,  the  spatial-spectral  scene,  and  the  atmospheric  seeing  parameter  for  the 
first  time.  The  algorithm  performance  was  measured  using  a  simple  simulation  and 
was  found  to  estimate  the  seeing  parameter  well.  In  the  future,  this  blind  estimator 
could  be  used  to  take  these  sensors  from  a  laboratory  environment  and  make  use  of 
existing  optical  telescopes  with  this  aberration  present.  Future  work  should  focus 
on  using  this  sensor  in  the  presence  of  atmosphere  and  an  experimental  test  of  this 
algorithm.  Further  work  could  also  be  done  on  designing  a  model  that  incorporates 
the  lower-bound  as  a  test  parameter  and  determining  the  limits  of  its  applicability 
to  broad  sensor  system  design  by  changing  the  parameters  and  testing  them  against 
actual  sensors. 
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clutter-to-noise  ratio,  see  CNR 
IF,  see  frequency 

independent  and  identically  distributed  data, 
see  i.i.d.  data 

jammer-to- noise  ratio,  see  JNR 

probability  of  false  alarm,  see  detection  prob¬ 
ability,  false  alarm  probability 
pulse  repetition  frequency,  see  PRF 
pulse  repetition  interval,  see  PRI 

radar  coordinate  system,  see  coordinate  sys¬ 
tem 

radar  cross  section,  see  RCS 

signal-to-interference  plus  noise  ratio,  see 
SINR 
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